AFG1.-TK-88-0303 


mCJILE-COET 


ID 

O  MODEL  DESCRIPTION  FOR  THE  SOCRATES  CONTAMINATION  CODE 

C\] 


Jamej  B. 


Elgin  and  Robert  L.  Sundberg 


Spectral  Sciences,  Inc. 
IV  South  Bedford  Street 
Burlington.  MA  01 803 


21  October  l! 38 


DTiC 


ZLECTE 
FEB  1  5  1983 

.D* 


Fine  Report 

23  S  ptember  !985  -  23  September  1988 


Approved  for  public  release;  distribution  unlimited 


AIR  FORCE  GEOPHYSI:  3  LABORATORY 
AIR  FORCE  SYSTEMS  JMMANI) 

UNITE 0  STATES  AIR  FORCE 

HANSCOM  AFB .  MASSACHUSETTS  01731-500 


« 


vs^Tif'anOKT 


REPORT  DOCUMENTATION  PaoE 

.ECU-  CLASSIFICATION 

lb.  RESTRICTIVE  MARKINGS 

-NONE  . . 

.  .  M'CUH  it  CuASSlflCATION  authority 

3.  DISTRIBUTION  /AVAILABILITY  OF  REPORT 

:o  DiC.ASMFiCAT'GN  /  DOWNGRADING  SCHEOULt 

approved  for  public  ielua.-.e; 
distil  bu  t  i  on  nn  1  i  in  i  t  cd 

PERFORMING  ORGANIZATION  REPORT  NUMBER(S) 


6a  NAME  Of  PERFORMING  ORGANIZATION 

Spectral  Sciences,  Inc. 


6c  ADDRESS  (C  ry,  State,  and  /IP  Code) 

111  South  Bedford  Street 
Burlington,  MA  01803-5128 


6b  OFFICE  SYMBOL 
(If  applicable) 

N/A 


5  MONITORING  ORGANIZATION  REPORT  NUMBER(S) 

AFCL-TR-88-  0303 


7a.  NAME  OF  MONITORING  ORGANIZATION 

A-ir  Force  Ceophysics  Laboratory 


7b.  ADDRESS  (City,  State,  and  /IP Code) 

Hanscom  Air  Force  Base 
Bedford,  MA  01731-5000 


8a.  NAME  OF  FUNDING /SPONSORING 
ORGANIZATION 


Electronic  Systems  Division 


8b.  OFFICE  SYMBOL 
(If  applicable) 

N/A 


8c  ADDRESS  (City,  State,  and  /IP Code) 

Air  Force  Systems  Command,  USAF 
Hanscom  AFB ,  MA  01731 


1 1  TITLE  (Include  Security  Classification) 


9.  PROCUREMENT  INSTRUMENT  IDENTIFICATION  NUMBER 

F19628-85-C-0195 


10.  SOURCE  OF  FUNDING  NUMBERS 


I  PROGRAM 
■  ELEMENT  NO. 

62 101F 


PROJECT 

NO. 


Model  Description  for  the  SOCRATES  Contamination  Code 


12  PERSONAL  AUTHOR(S) 

James  li.  Elgin  and  Robert  L.  Sundber 


13a  TYPE  OF  REPORT  ll 3b.  TIME  COVERED  114.  DATE  OF  REPORT  (Year.  Month.  Day)  (l  5.  PAGE  COUNT 

Final  Report  I  FROM85SEE>23  T088SFP23  1  21  October  1988  I  128 


16  SUPPLEMENTARY  NOTATION 


17 

COSATI  CODES 

FIELD 

|  GROUP  1 

1  SUB-GROUP 

18.  SUBJECT^T^RMS  (Continue  on  reverse  if  necessary  and  identify  by  block  number) 

Contamination  Monte  Carlo  Space  Shuttle,  {r 

Rarefied  Flows 


19  ABSTRACT  (Contlnde  on  reverse  if  necessary  and  identify  by  block  number) 


Ti\h  SOCRATES  contamination  model  is  described  at  length.  The  model 
aiiuws  Tor  unsteady  or  steady  simulation  or  contamination  aboard  the  space 
snuitie  orbiter  via  the  direct  simulation  Monte  Carlo  method.  The  basis 
tor  trie  model  is  discussed,  arid  sample  calculations  are  given  for  an  RCS 
engine  tiring  at  200,  250,  ano  300  kilometer  attitudes.  The  dependence  of 
renirtt  flux  of  scattered  species  on  exhaust  species  molecular  weight  is 
e :  in:  Miatea. 

X 


V 


20  DISTRIBUTION /AVAIIA8IUTY  OF  ABSTRACT 

□  UNCLASSIFIED/UNLIMITED  □  SAME  AS  RPT.  □  OTIC  USERS 


22a  NAME  OF  RESPONSIBLE  INOIVIOUAL 

Dr.  Shu  Lai 


OO  FORM  1 473, 84  MAR  83  APR  adition  may  ba  used  until  axhaustad 

All  ottiar  edition*  are  obtolata 


121.  ABSTRACT  SECURITY  CLASSIFICATION 
□  OTIC  USERS  I  UNCLASSIFIED _ 


22b.  TELEPHONE  (Include  Area  Code)  1 22c.  OFFICE  SYMBOL 

(617)  377-2 


IllilO 


UNCLASSIFIED 


CONTENTS 


1  OVERVIEW  OF  THE  SOCRATES  COOK .  1 

2  GAS  MODEL  AND  EQUILIBRIUM  PROPERTIES  .  4 

2.1  Preliminary  Equilibrium  Gas  Relations  .  4 

2.2  Analytical  Form  oT  the  Collision  Cross  Section  .  6 

2.3  Equilibrium  Reference  Properties  for  a  Multi- 

Component  Gas .  8 

2.4  Internal  Energy  Model  .  10 

3  INTERNAL  REPRESENTATION  .  11 

3.1  State  Vector .  11 

3.2  Reduction  to  a  Reasonable  Number  of  Simulated 

Molecules .  11 

3.3  Internal  Scales . 12 

3.4  Weigricmg  Factors  .  33 

4  GRID  COORDINATES  AND  GRID  STRUCTURE . 14 

5  COLLISION  MECHANICS  .  16 

5.1  Relations  for  Elastic  Collisions  .  16 

5.2  Effect  ol  coordinate  System .  .  .  18 

5.3  Simulation  of  Inelastic  Collisions .  1!) 

5.4  Collisions  Between  Molecules  with  Distinct 

Weighting  Factors . 21 

5.5  Reactive  Collisions  .  21 

5.5.1  Types  of  Reactive  Collisions  .  22 

5.5.2  Reactive  Collision  Probability  .  23 

5.5.3  Options  for  Simulating  Reactive 

Collisions . 24 

6  COLLISION  SAMPLING  IN  A  MULTI -COMPONENT  VHS  GAS  ...  26 

6.1  General  Considerations  and  Approach  .  26 

6.2  Collision  Sampling  for  a  Single  Component  Gas  27 

6.2.1  Collision  Pair  Selection . 27 

6.2.2  Collision  Time  Counter  for  a  Single 

Component  Gas . 28 

6.3  Collision  Class  Sampling  in  Gas  Mixtures  ....  29 

6.4  Global  Collision  Sampling  in  a  Gas  Mixture  ...  30 

6.4.1  Global  Collision  Time  Counter  .  30 

6.4.2  Collision  Pair  Selection  in  Multi¬ 
component  Mixtures . 32 


i  i  i 


6.4.3  Summary  of  Collision  Sampling  in  Multi¬ 
component  Mixtures . 32 

6.5  Deviations  from  the  General  Procedure  .  33 

6.5.1  Cell  Specific  Atm . 33 

6.5.2  Relaxation  of  Q,liax . 34 

6.5.3  Maximum  Time  Counter . 34 

6.5.4  Separation  of  Major  ana  Minor  Species  .  .  35 

PROCEDURES  FOR  COLLISION  DOMINATED  FLOW  .  36 

7.1  Collision  Cutoff  Approach  .  36 

7.2  Equilibrium  Aftermath  Approach  .  37 

7.2.1  Conserved  Quantities  .  38 

7.2.2  Center  of-Mass  Velocity  Distribution  .  .  39 

7.2.3  Molecular  Relative  Velocity 

Distribution  .  41 

7.2.4  Translational  Energy  of  Relative  Motion  .  41 

7.2.5  Determination  of  Temperature  .  42 

7.3  The  Number  of  Collisions  Requirea  to  Achieve 

Equilibrium . 43 

7.4  Method  Comparison  .  44 

MOLECULAR  TRANSLATIONS  .  44 

8.1  Molecular  Translations  in  Cartesian 

Coordinates . 44 

8.2  Molecular  Cloning  .  45 

INITIAL  AND  BOUNDARY  CONDITIONS  .  47 

9.1  Initial  Conditions  .  47 

9.1.1  Number  of  Simulated  Molecules  and  Their 

Weighting  Factors  .  47 

9.1.2  Initial  Positions  .  48 

9.1.3  Initial  Velocity  Components  .  48 

9.1.4  Initial  Internal  Energies  .  49 

9.2  Source  Boundary  Conditions  .  49 

9.3  Atmospheric  Boundary  Condition  .  51 

9.3.1  Molecular  Flux  Across  a  Surface  Element  .  51 

9.3.2  Incoming  Molecular  Velocity  Components  52 

STATISTICAL  SAMPLING  OF  OUTPUT  .  53 

10.1  General  Considerations  .  53 

10.2  Sampling  of  Instantaneous  Volumetric  Output 

Quantities . 54 

10.3  Sampling  of  Time  Averaged  Output  Quantities  .  .  57 

SURFACE  DEFINITIONS  AND  INTERACTIONS  .  57 

11.1  Shuttle  Representation  .  58 


11.2  Determination  of  Surface  Intersections  ....  62 

11.3  Surface  Reflections  .  66 

11.4  Surface  Statistics  .  66 

11.5  Interface  of  Shuttle  Model  with  Calculationai 

Grid . 66 

;1.6  Back-to-Back  Surfaces  .  67 

12  SAMPLE  CALCULATIONS  .  67 

12.1  Case  Descriptions . 67 

12.2  Contamination  Cloud  Results  at  200  Kilometers  .  70 

12.3  Surface  Contamination  at  200  Kilometers  ....  86 

12.4  Results  at  250  anu  300  Kilometers . 91 

13  CONCLUSIONS . 99 

14  REFERENCES . 99 

APPEND f  X  A . A-l 


v 


- 

ILLUSTRATIONS 

i 

A  Schematic  Representation  of  the  Major  Elements  of 
Shuttle  Contamination  Problem  . 

3 

- 

- 

A  Diagram  of  the  Basic  Solution  Proceoure  Utilized 
for  Steady  State  Solutions  in  the  SOCRATES 
Contamination  .Monel . 

4 

3 

A  Schematic  Showing  the  Basic  Design  of  the  Socrates 
Cell  Structure  Which  Uses  Cartesian  Coordinates  with 
Uneven  Spacing  . 

16 

4 

A  Frontal  View  of  the  Crude  Shuttle  Model  Designed 
for  Testing  of  the  SOCRATES  Monel  . 

59 

5 

A  Top  View  of  the  Crude  Shuttle  Model  Designed  for 
Testing  of  the  SOCRATES  Model  . 

60 

- 

6 

A  Side  View  of  the  Crude  Shuttle  Model  Designed  for 
Testing  of  the  SOCRATF.S  Model . 

61 

- 

7 

An  Illustration  of  the  Quantities  Used  to  Calculate 
a  Molecular  Intersection  with  a  Rectangular  Plate  .  . 

63 

- 

f! 

A  Schematic  of  the  Sample  Shuttle  Problem  Showing 
the  Coordinate  System  and  Orientation  of  the 
Calculation  . 

66 

9 

A  Contour  Plot  of  the  Scattered  H2  Number  Density 
at  an  X  Location  of  0  Meters  For  the  200  Kilometer 

Case . 

72 

!  0 

A  Contour  Plot  of  the  Scattereu  H2  Number  Density 
at  an  X  Location  of  500  Meters  For  the  200  Kilometer 

73 

1 1 

A  Contour  Plot  of  the  Scattered  H2  Number  Density 
at  an  X  Location  of  1500  Meters  For  the  200  Kilometer 
Case . 

74 

12 

A  Contour  Plot  of  the  Scattered  H2  Number  Density 
at  a  2  Location  of  1000  Meters  For  the  200  Kilometer 
Case . 

75 

- 

vi 

’.3  A  Contour  Plot  of  the  Scattered  C02  Number  Density 

at  a  Z  Location  ot  1000  Meters  For  trie  200  Kilometer 
Case . 76 

:4  The  Upstream  Density  Decay  for  the  Three  Species  at 
a  L  Location  of  ]  Kilometer  Above  the  Shuttle  for 
tiie  200  Kilometer  Case . 77 

15  A  Contour  Plot  of  the  Z  Velocity  Component  for  the 
Scattered  H2  at.  an  X  Location  of  0  Meters  For  the 

200  Kilometer  Case . 78 

16  A  Contour  Plot  of  the  Translational  Temperature  for 
the  Scattered  H2  at.  an  X  Location  of  0  Meters  For 

the  200  Kilometer  Case . 79 

17  A  Contour  Plot  of  the  Z  Velocity  Component  for  the 
Scattered  C02  at.  an  X  Location  of  0  Meters  For  the 

200  Kilometer  Case . 81 

1 8  A  Contour  Plot  of  the  Translational  Temperature  for 
trie  Scattered  C02  at  an  X  Location  ol  0  Meters  For 

the  200  Ki  lometer  Case . 82 

19  Trie  Effect  of  Exhaust  Species  Molecular  Weight  on 

the  Normalized  Scattered  Species  Density  in  the 
Vicinity  of  the  Shuttle . 83 

20  The  Effect  of  Exhaust  Species  Molecular  Weight  on 

the  Scattered  Species  Z  Velocity  Component  in  the 
Vicinity  of  tiie  Shuttle . 84 

2:  Tiie  Effect  of  Exhaust  Species  Molecular  Weight  on 

trie  Scattered  Species  Translational  Temperature  In 
trie  Vicinity  ot  the  Shuttle . 85 

22  The  Return  Flux  of  Scattered  Hg  in  tiie  Y-Z  Plane 
as  a  Function  of  Azimuthal  Anglo,  in  tire  Vicinity 

of  tho  Shuttle,  for  the  200  Km  Case . 88 

23  The  Return  Flux  of  Scattered  H20  in  the  Y-Z  Plane 
as  a  Function  of  Azimuthal  Angle,  in  the  Vicinity 

of  the  Shuttle,  for  tiie  200  Km  Case . 89 


vii 


24 


The  Return  Flux  ot  Scattered  C02  in  the  Y-Z  Plane 
us  a  Function  of  Azimuthal  Angle,  in  the  Vicinity 
of  the  Shuttle,  for  the  200  Km  Case . 90 

25  The  Return  Flux  of  Scattered  H2  in  the  Y-Z  Plane 
as  a  Function  of  Azimuthal  Angle,  in  the  Vicinity 

of  the  Shuttle,  for  the  250  Km  Case . 92 

26  The  Return  Flux  of  Scattered  H20  in  the  Y-Z  PJane 
as  a  Function  of  Azimuthal  Angle,  in  the  Vicinity 

of  the  Shuttle,  for  the  250  Km  Case . 93 

27  The  Return  Flux  of  Scattered  C02  in  the  Y-Z  Plane 
as  a  Function  of  Azimuthal  Ancle,  in  the  Vicinity 

of  the  Shuttle,  for  the  250  Km  Case . 94 

28  The  Return  Flux  of  Scattered  H2  in  the  Y-Z  Plane 
as  a  Function  of  Azimuthal  Angle,  in  the  Vicinity 

of  the  Shuttle,  for  the  300  Km  Case . 95 

29  The  Return  Fiux  of  Scattered  H20  in  the  Y-Z  PJane 
as  a  Function  of  Azimuthal  Angle,  in  the  Vicinity 

of  the  Shuttle,  for  the  300  Km  Case . 96 

30  The  Return  Flux  of  Scattered  C02  in  the  Y-Z  Plane 
as  a  Function  of  Azimuthal  Angle,  in  the  Vicinity 

of  the  Shuttle,  for  the  300  Km  Case . 97 

31  The  Dependence  of  Maximum  Return  Flux  of  the 
Atmospherically  Scattered  MoJecules  on  Ambient 

Number  Density  for  the  200,  250,  and  300  Km  Runs  .  .  98 

A-J  A  Plot  of  the  Chi-Square  Probability  Density  Function 

for  V  Equal  to  1,  2,  and  3 .  A-2 

A-2  A  Representation  of  the  Transformed  Chi-Square 

Distribution,  p(Z),  for  v  =  50  (Solid  Line)  ....  A-10 

A-3  A  Representation  of  the  Transformed  Chi-Square 

Distribution,  h(W),  for  v  -  1  (Solid  Line)  .  A-12 


A-4 


A  Comparison  of  the  Exact  and  Approximate  Acceptance 
Probabilities  for  v-1  . 


A- 16 


A -5  A  Comparison  of  the  Exact  Chi-Square  Distribution, 
t ( X ; v )  with  the  Approximate  Distribution. 

t  a  ( X ;  V ) .  A-17 

TABLES 

1  Scaling  Factors  Used  for  the  Internal  Representation 

of  Quantities  in  the  Socrates  Code . 13 

2  Kxnaust  Composition  for  Sample  Calculations  .  69 

3  Atmospneric  anu  Caicu.ai.i oriai  Parameters  for  Sample 

Calculations  .  70 


ix 


OVERVIEW  OF  THE  SOCRATES  CODE 


1  . 


The  airect  simulation  Monte  Carlo  method,  as  pioneered  by  G.  A.  Bird,1 
provides  a  powerful  technique  for  the  simulation  of  real  gas  flows.  It 
bridges  the  gap  between  continuum  and  free  molecular  flow,  retaining 
validity  in  either  extreme.  It  can  be  used  to  describe  complex  mixtures, 
including  effects  of  chemical  reactions,  heat  conduction,  viscosity  and 
diffusion  for  flows  in  three  dimensions.  This  report  describes  the 
application  of  this  technique  to  the  contamination  problem,  considering 
flow  fields  created  by  the  interaction  of  a  spacecraft  with  the  atmosphere. 
The  resultant  model  has  been  named  the  SOCRATES  code,  which  is  an  acronym 
for  Shuttle  Orbiter  Contamination  Representation  Accounting  for  Transiently 
Emitted  Species. 

Contamination  of  instruments  on  the  space  shuttle  orbiter  is  an  issue 
of  major  concern.  The  shuttle  gives  off  matter  through  surface  outgassing, 
via  various  thrusters,  and  from  flash  evaporators.  At  altitudes  where  the 
atmospheric  mean  free  path  is  comparable  to  or  less  than  shuttle 
dimensions,  the  deposition  back  onto  shuttle-borne  instruments  will  be 
largely  determined  by  the  multiple  collision  environment  surrounding  the 
shuttle.  Even  at  higher  altitudes,  this  may  be  the  dominant  source  of 
contaminants  for  some  portions  of  the  shuttle.  In  addition  to  physical 
contamination  of  shuttle  surfaces,  "radiation  contamination"  is  also  a 
potential  problem  as  gases  surrounding  the  shuttle  collide  at  high  speed 
with  atmospheric  molecules.  These  energetic  collisions  can  lead  to 
vibrational  excitation  and  subsequent  radiative  decay.  A  similar  issue  of 
some  concern  is  the  presence  of  ions  in  the  vicinity  of  the  shuttle  which 
can  (possibly)  be  produced  via  the  critical  ionization  velocity  effect. 
Ions  in  the  shuttle  environment  may  remain  there  for  some  time  due  to 
eiectric  field  forces,  and  radiative  recombination  is  another  potential 
source  of  radiation  contamination.  The  situation  is  depicted  schematically 
in  Figure  1 . 

Spectral  Sciences,  Iric.  (SSI)  developed  an  initial  version  of  this 
three-dimensional  Monte  Carlo  model  of  the  flow  field  about  the  shuttle  so 
that  the  contamination  can  be  accurately  characterized  and  understood.  A 
comprehensive  model  of  the  contaminant  field  surrounding  the  space  shuttle 
orbiter  is  crucial  to  the  design  of  experiments  which  are  to  fly  on  the 


1.  Bird,  G.  A.,  Molecular  Gas  Dynamics,  Clarendon  Press,  Oxford  (1976). 


shuttle  and  to  the  development  of  procedures  for  minimizing  the 
contamination .  The  code  is  designed  in  a  highly  modularized  fashion  so 
that  additional  physical  and  geometric  complexity  can  be  added  as  deemed 
necessary  without  requiring  major  rewriting  of  the  model. 

The  basic  caiculational  technique  is  well  described  by  its  originator 
in  Ret.  1.  However,  there  have  been  significant  extensions  of  the  method 
since  the  publication  of  Ref.  1.  The  present  purpose  is  to  describe  how 
the  technique  is  implemented  in  SOCRATES;  but  elementary  concepts  and 
relations  which  are  essential  to  a  coherent  explanation  are  included  here 
also. 


The  direct  simulation  Monte  Carlo  method  involves  storing  a  discrete 
number  of  molecules  (via  their  velocities,  positions,  and  other  pertinent 
information)  in  a  computer.  The  solution  region  is  broken  up  into  a  number 
of  separate  cells,  and  the  solution  is  stepped  forward  in  time  in  a  two- 
stage  process.  First,  the  molecules  are  advanced  along  their  trajectories 
by  art  amount  appropriate  to  their  velocity  ana  a  time  increment,  Atm.  In 
this  first  stage  some  molecules  will  leave  the  solution  region,  and  some 
will  be  introduced  as  determined  by  the  boundary  conaitions  for  a 
particular  problem.  The  second  stage  is  to  simulate  collisions  in  each 
cell  appropriate  to  Atm  so  that  collision  frequencies  are  properly 
simulated.  A  basic  hypothesis  of  the  method  is  that  if  the  time  step  is 
made  small  enough,  the  processes  of  translations  and  collisions  can  be 
uncoupled  in  this  manner. 

Periodically,  the  solution  is  sampled  by  accumulating  statistical  sums 
of  number  densities,  velocities  and  other  basic  properties.  The  solution 
is  run  repeatedly  until  statistical  deviations  are  reduced  to  a  desired 
limit,  and  then  physically  meaningful  output  quantities  are  computed  from 
the  statistical  sums.  The  number  of  molecules  represented  is  typically 
many  thousand  at  a  time,  which  is  vastly  fewer  than  the  number  occurring  in 
virtually  all  real  flows.  Hence,  the  construction  of  a  dynamically  similar 
flow  to  be  simulated  in  the  computer  is  an  essential  feature  of  the  method. 

The  logical  flow  of  the  solution  procedure  is  shown  in  Figure  2,  which 
includes  the  steps  described  above.  The  following  sections  describe  in 
detail  the  Implementation  of  each  of  the  boxes  shown  in  Figure  2  and  the 
application  of  the  code  to  some  sample  problems. 
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Figure  1.  A  Schematic  Representation  oi  the  Major  Elements  of  the 
Shuttle  Contamination  Problem. 


SAMPLE  CURRENT  MOLECULAR  PARAMETERS 


ND 


GENERATE  OUTPUT  FROM  ACCUMULATED  STATISTICS 


Figure  2.  A  Diugram  of  the  Basic  Solution  Procedure  Utilized  lor 
Steady  State  Solutions  in  the  SOCRATHS  Contamination 
Model . 

2.  GAS  MODEL  AND  EQUILIBRIUM  PROPERTIES 

2.1  Preliminary  Equilibrium  Gas  Relations 

The  far  field  equilibrium  state  has  properties  which  are  of  relevance 
to  the  flow  field  interaction  problem  to  be  solved.  Length  and  velocity 
scaies  for  the  problem  are  obtained  from  the  far  field  and  used  to 
non-aimensional ize  the  internal  code  variables.  Even  if  this  were  not 
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aone,  it  provide  an  important  comparison 

velocities,  collision  frequencies,  etc. 


case  for  densities, 


?or  a  rest  gas  in  equilibrium,  the  normalized  distribution  function 
ror  trie  relative  speed,  cr,  between  molecules  of  species  i  and  molecules  of 


species  j  is  given  by 


2 


*  i  1  < c  r ) 


,2  3/2  , 


(4craij  exp(-aijC“) 


(i: 


where 


‘ij 


i.l 


2*0^ 


(2] 


ano 


M, j  is  the  reduced  mass  of  the  pair;  i.e., 


mim  j 

mj+mj 


(3) 


with  id,  ano  mj  representing  the  masses  of  the  two  species.  In  these 
relations,  Ta  is  the  far  field  temperature  and  RQ  is  the  universal  gas 
constant.  {Rq  is  used  instead  of  Boltzmann's  constant  since  the  molecular 
masses  will  be  consistently  represented  in  atomic  mass  units  rather  than 
grams.)  The  available  translational  coilisional  energy  between  the  two 
molecules,  Ec,  is  given  by 


(4) 


2.  Chapman,  S.  and  Cowling,  T.  G.  ,  The  Mathematical  Theory  of  Non- 
Uniform  Gases,  3rd  ed.,  Cambridge  University  Press,  Cambridge,  86 
( 1070)  . 


2.2  Analytical  Form  of  the  Collision  Cross  Section 


Whenever  the  direct  simulation  Monte  Carlo  method  is  applied,  it  is 
necessary  to  make  trade  oils  between  accuracy  and  simplicity  in  molecular 
models.  It  does  no  good  to  use  a  complex  molecular  potential  surface  and 
then  find  that  reasonable  computer  run  times  result  in  very  large 
statistical  fluctuations  in  the  output.  Since  the  final  output  will 
reflect  errors  in  the  statistics  as  well  as  errors  in  the  models,  there  is 
a  strong  impetus  to  use  models  which  contain  the  essential  physics,  but 
which  can  be  applied  in  a  computationally  efficient  manner.  The  current 
state-of-the-art  is  the  Variable-Hard-Sphere  (VHS)  model  3  In  this  model 
mojecules  have  a  collision  cross  section  which  varies  as  an  inverse  power 
of  the  available  collision  energy.  Hence,  if  cr^j  is  the  collision  cross 
section  for  collisions  of  species  i  with  species  j,  then  cr^  is  given  by  a 
relation  of  the  form 

ffij  =  AuC  •  <5> 

where  Axj  is  a  constant  coefficient.  It  follows  that  the  effective 
diameter  for  molecules  of  species  i,  dj  ,  is  implicitly  defined  as  a 
function  of  available  collision  energy  by  the  relation 

CU  =  iril'f  =  AUE^  (6) 
Ai:t  can  be  determined  from  a  reference  cross  section  ano  velocity  via 

Aii  =  •-cri  x  (micr/4  )uJref  •  ^ 

If  a  reference  cross  section  is  given  for  a  reference  temperature  rather 
than  a  reference  velocity,  then  the  usual  choice  for  the  reference  velocity 
is  that  velocity  which  has  a  collision  energy  equal  to  the  mean  collision 
energy  occurring  in  collisions  at  the  reference  temperature. 
Mathematically,  this  is  equivalent  to 


3.  Bird,  G.  A.,  "Morite-Cario  Simulation  in  an  Engineering  Context”, 
Proceedings  of  the  12th  International  Symposium  on  Rarefied  Gas 
Dynamics,  74,  Progress  in  Astronautics  and  Aeronautics,  AIAA ,  New 
York  (1981 )  .  .  ~  ~  ~ 
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(8) 


2  <cr3<rii> 

(cr  ref  =  "<^nT 

where  the  angle  brackets  indicate  averages  taken  over  the  distribution 
function  given  in  Eq.  (1)  evaluated  for  and  T^Tpg^.  Equation  (8) 

can  be  simplified  to  give 


(cr^ref 


4(2  -  t>)RoTref 


(0) 


For  simulations  involving  a  large  number  of  species,  reference  cross 
sections  are  frequently  not  available  for  all  possible  collision  pairs.  In 
this  case  it  is  possible  to  specify  for  self-collisions  only,  and  then 
use  Eq.  (6)  to  get  a  molecular  diameter  as  a  function  of  collision  energy. 
Tfien.  applying  the  relation 

ffij  =  ff[(dj  +  a3)/2\2  .  (10) 
the  coefficient  in  Eq.  (5)  for  interspecie  collisions  is  given  by 


Ai j  =  i-|(^Aii  +  ^Ajj  H2 


(ID 


for  the  internal  workings  of  a  Monte  Carlo  code,  it  is  usually  more 
convenient  to  express  the  collision  cross  section  as  a  function  of  the 
relative  collision  velocity  rather  than  the  collision  energy.  This  is 
simply  achieved  via  the  relation 


U 


„  b  c~2u 
BiJcr 


where 


3ij  ~  Aij(uij/2)  “ 


(12) 


(13) 


The  parameter  u  can  be  related  to  t),  the  exponent  of  distance  in  an 

O 

inverse  power  interraolecular  force  law  via  the  relation 


(0 


(14) 


-  7  - 


Hence,  hard  sphere  molecules  (for  which  T)  goes  to  infinity)  are  represented 
by  a  equal  to  zero.  There  is  a  substantial  body  of  evidence,  however,  that 
the  effective  size  of  molecules  does  indeed  decrease  with  increasing 
collision  energy,  so  a  positive  value  of  ci  is  usually  a  belter  choice,  o 
can  be  determined  from  molecular  beam  data,  or  from  its  macroscopic 
implications.  For  example,  if  s  is  the  exponent,  for  the  variation  of  the 
viscosity  coefficient  with  temperature,  then  it  can  be  shown3  that 

s  -  «  +  0 . 5  (15) 

so  a  measurement  of  the  temperature  dependence  of  the  viscosity  coefficient 
serves  as  an  indirect  determination  of  o. 

In  oraer  to  incorporate  the  model  for  internal  energy  transfer  to  be 
discussed  in  Section  5,  it  is  necessary  that  u  be  assumed  the  same  for  all 
interactions.  This  represents  one  of  the  major  restrictions  in  the  current 
state  or  monel ing. 

Although  the  sizes  of  molecules  are  allowed  to  vary  in  the  VHS  model 
m  deciding  whether  or  not  a  collision  is  to  occur,  when  a  collision  does 
occur  the  post  collision  velocity  components  are  computed  as  if  it  were  a 
hard  sphere  collision  (see  Section  5).  This  results  in  a  substantial 
computational  simplification  and  yet  retains  good  agreement  with  the 
macroscopic  predictions  of  the  more  exact  model.3  (See  Ref.  1  for  a 
niscussion  of  molecular  scattering  for  general  power  law  potentials.) 

2,3  Equilibrium  Reference  Properties  for  a  Multi-Component 
Gas 

One  advantage  of  the  VHS  model  is  that  the  molecules  have  a  well 
defined  cross  section,  so  it  is  possible  to  define  a  mean  free  path  without 
putting  limitations  on  the  minimum  deflection  angle  that  is  considered.  As 
is  the  general  case  for  multi-component  gases,  however,  each  component  has 
its  own  mean  free  path,  and  the  overall  mean  free  path  for  the  mixture  must 
be  defined  as  a  weighted  average  of  the  mean  free  paths  of  the  individual 
species.  The  somewhat  cumbersome  relations  required  to  calculate  the 
overall  mean  free  path  are  given  here.  It  should  be  noted  that  the  mean 
free  path  is  calculated  only  once  lor  a  given  problem,  so  the  computational 
effort  required  to  evaluate  it  is  completely  negligible. 
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An  individual  molecule  of  species  i  will  suffer  collisions  with 
molecules  of  species  j  with  a  frequency  vj  given  by 


vi  =  njto<ffi  jcr> 


where  nJaj  is  the  number  density  of  species  j  and  «r^jCr>  is  the  average 
product  of  cross  section  times  relative  velocity  for  the  two  species, 
oocairted  by  integrating  over  the  distribution  function  given  in  Eq.  (1). 
When  this  operation  is  performed,  the  result  is 


vl  =  2BijnjcDr<2  - 


wnere  T  denotes  the  gamma  function. 

The  total  collision  frequency  for  an  individual  molecule  of  species  i, 
Vj .  is  obtained  by  summing  Eq.  (16)  over  all  species,  i.e., 


;i  =  IV1 


ana  the  mean  free  path,  X j ,  for  molecules  of  species  i  is  given  by 


-•VSKoT^/firmj) 


wnere  <Cj>  is  the  mean  molecular  speed  for  species  i  molecules.  The  mean 
free  path  for  the  gas  mixture,  X^,  is  then  defined  as  the  number  density 
weigntea  average  of  the  X^  via 


9  . 

r*  n<_X, 


V  ni®  i 
=  >  — 


where  nffl  is  the  total  number  density: 


r 
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A  useful  velocity  scale  is  given  by  vg,  defined  by 


vs  =  -V2R0TCB/<m>  . 

where  <m>  is  the  reference  mean  molecular  weight,  i.e.. 


(22) 


<m> 


P 

V" 

; 


nicon'i 


i=l 


(23) 


vs  is  the  most  probable  molecular  speed  for  molecules  of  the  mean  molecular 
weight  at  the  reference  temperature. 


2,4  Internal  Energy  Model 


The  current  state  of  modeling  for  internal  energy  effects  in  Monte 
Carlo  flow  field  simulations  is  the  phenomenological  model  of  Borgnakke  and 
Larsen.4  In  this  model,  transfer  of  energy  between  internal  and 
translational  modes  is  allowed,  but  it  is  necessary  to  assume  that  each 
species  has  a  fixed  number  of  internal  degrees  of  freedom,  .  This  is 
equivalent  to  assuming  a  constant  specific  heat,  Cp ^ ,  for  each  species 
which  can  be  related  to  the  number  of  internal  degrees  of  freedom  via 


=  2 


Jpi  i 


R 


-  5 


0 


(24) 


Alternatively,  Cj  can  be  related  to  the  ratio  of  specific  heats  for  species 
i,  Tj,  by  the  relation 


The  interchange  of  internal  and  translational  energy  will  be  discussed  in 
Section  5,  and  the  selection  of  initial  conditions  will  be  discussed  in 
Section  9. 

4.  Borgnakke,  C.  and  Larsen,  P.  S. ,  "Statisticai  Collision  Model  for 
Monte  Carlo  Simulation  of  Polyatomic  Gas  Mixture",  Journal  of 
Computational  Physics,  18,  405  (1975). 


3 .  INTERNAL  REPRESENTATION 

3 . 1  State  Vector 

Each  simulated  molecule  in  the  SOCRATES  code  is  represented  by  a  state 
vector  which  comprises  all  of  the  information  the  code  has  with  regard  to 
that  particular  molecule.  The  state  vector  has: 

•  Position  elements  defining  the  location  of  the  molecule  in 
Cartesian  coordinates. 

■  Three  velocity  elements,  giving  the  corresponding  velocity 
components  in  the  same  coordinate  system. 

•  A  value  for  the  internal  energy  (usually  rotational)  of  the 
molecule.  Note  that  the  basic  model  does  not  discriminate 
between  internal  modes  for  a  particular  species.  This  can  be 
done,  if  desired,  by  introducing  separate  species  for  the 
distinct  modes. 

•  An  indicator  identifying  the  molecular  species.  This 
indicator  in  turn  implies  all  of  the  properties  associated 
with  that  species  (molecular  weight,  number  of  internal 
degrees  of  freedom,  name,  etc.). 

•  An  indicator  giving  the  computation  cell  in  which  the  molecule 
currently  resides.  (This  could  be  calculated  from  its 
position,  but  it  is  needed  so  often  in  the  calculation  that 
the  extra  storage  location  is  justified  by  the  increase  in 
efficiency. ) 

•  A  time  element,  giving  the  time  at  which  the  molecule  will 
strike  a  solution  surface  element  if  it  continues  on  its 
current  trajectory.  (See  the  discussion  in  Section  11.) 

3.2  Reduction  to  a  Reasonable  Number  of  Simulated  Molecules 

It  is  clearly  impossible  to  run  a  computer  simulation  with  anywhere 
near  the  same  number  of  molecules  that  exist  In  the  actual  flow  problem. 
The  adjustment  that  is  maae  to  make  the  simulation  possible  is  to 
artificially  increase  the  cross  section,  and  decrease  the  number  density, 
by  the  same  large  factor.  ft  is  the  product  of  number  density  and  cross 
section  which  determines  the  collision  frequency  for  a  given  molecule,  and 
it  is  the  collision  frequency  which  must  be  correctly  simulated  if  the 
correspondence  between  the  real  and  simulated  flows  is  to  be  accurate. 
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This  is  an  essential  feature  ol  the  direct  simulation  method  which  has  not 
always  been  adequately  emphasized.  It  means  that  the  internal  scaling 
factors  do  not  proceed  on  a  strictly  dimensional  basis.  For  example,  the 
scaling  factor  for  cross  sections  is  not  the  square  of  the  scaling  factor 
for  lengths. 

3.3  Internal  Scales 

Many  problems  are  more  reasonably  handled  if  the  internal  calculations 
are  carried  out  with  scaled  or  dimensionless  values.  This  avoids  possible 
problems  such  as  numerical  overflow  which  can  cause  an  execution  time 
error.  Such  errors  can  be  particularly  insidious  and  difficult  to  locate 
in  a  code  whose  very  essence  involves  the  random  combination  of  numbers. 

The  output  is  produced  in  physically  meaningful  dimensional  form. 
Hence,  the  scaling  that  is  discussed  here  is  irrelevant  (or  nearly  so)  to 
the  interpretation  of  code  output;  it  is  strictly  a  matter  of  the  internal 
representation . 

The  choices  for  length  and  velocity  scales  are  X.^  and  vs  as  defined  in 
Section  2,  which  are  used  to  non-dimensionalize  the  position  and  velocity 
elements  of  the  state  vector.  There  is  no  need  to  provide  further 
non-dimensional ization  of  mass  beyond  representing  them  in  atomic  mass 
units,  so  none  is  provided.  Hence,  the  scaling  factor  for  energy  is  just 
v^ ,  which  is  used  to  non-dimensionai tze  the  internal  energy  element  of  the 
state  vector. 

Number  densities  are  scaled  with  respect  to  the  far  field  ambient 
number  aensity,  n^,.  which  leaves  only  the  cross  section  scaling  factor  to 
be  determined.  This  factor  follows  from  the  condition  of  flow  similarity, 
which  requires  that  the  probability  of  a  molecule  suffering  a  collision  in 
traveling  a  given  path  length  be  accurately  simulated.  This  dimensionless 
proDability  can  be  expressed  as  the  product  of  a  cross  section  times  a 
number  density  times  a  path  length  (at  least  for  small  enough  path 
lengths),  and  it  is  required  that  this  product  be  the  same  for  dimensional 
and  scaled  representations.  This  implies  that  the  product  of  the  scaling 
factors  for  these  three  quantities  be  unity  and,  therefore,  that  the  cross 
section  scaling  factor  be  l/fn^X.^).  The  internal  scaling  factors  used  in 
SOCRATES  are  summarized  in  Table  1. 
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Table  1. 


Scaling  Factors  Used  for  the  Internal 
Representation  of  Quantities  in  the  Socrates 
Coae.  All  Variables  are  Defined  in  Section 
2. 


PROPERTY 

SCALING  FACTOR 

Length  . 

Velocity  . 

..  vs 

Time  . 

..  \ra/vs 

Number  Density  .. 

■  •  n® 

Mass  . 

.  .  a .  m .  u . 

Energy  . 

. .  (a.m.u. )v^ 

Cross  Sect  ion  .  .  . 

3.4  Weighting  Factors 


Statistical  weigntir.g  factors  are  a  crucial  element  of  a  successful 
Monte  Carlo  simulation,  allowing  trace  species  to  be  described  with 
reasonable  accuracy.  The  weighting  factor  is  the  number  of  "real" 
molecules  that  correspond  to  each  "simulated"  molecule.  A  "simulated" 
moiecuie  corresponds  to  one  molecule's  worth  of  storage  (one  state  vector) 
allocated  in  the  program,  and  the  weighting  factor  is  its  statistical 
weight.  So,  for  example,  the  total  number  density  in  a  cell,  n(.eJl  can  be: 
expressed 


"cel  i 


y 


a 

V 


1=1 


(26) 


where  Ni  indicates  the  numDer  of  simulated  molecules  of  species  i  in  the 
ceil.  is  the  weighting  factor  for  that  species  in  that  cell,  V  Is  the 
ceil  voiume,  ar.o  p  is  the  number  of  species.  The  product  NXWX  that  appears 
m  Eq.  (26)  is  termed  the  number  of  "real"  molecules  of  species  i  in  the 
ce:i  .  Note  that  nfi  ^  as  calculated  by  Eq.  (26)  is  a  scaled  value;  it 
would  have  to  be  multiplied  by  n  ,  as  shown  in  Table  1,  to  become  a 
dimensional  evaluation  of  the  number  density. 


The  weighting  factors  used  in  SOCRATES  are  dependent  on  ceil  anti 
species.  Hence,  flow  fields  where  a  given  species  is  much  more  dominant  in 
one  portion  of  the  solution  region  than  another  can  be  accurately 
represented. 
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A  critical  error  that  can  occur  in  Monte  Carlo  codes  is  to  have  the 
number  or  simulated  molecules  exceed  the  dimensioned  limit  of  the  code.  On 
the  other  hand,  it  is  generally  desirable  to  have  as  many  molecules  as  is 
feasible  to  obtain  good  statistics.  Resolution  of  these  conflicting 
considerations  is  complicated  by  lack  of  a  priori  knowledge  of  what  the 
species  number  densities  will  be  as  a  function  of  space  and  time.  The  way 
the  resolution  is  achieved  is  by  a  dynamic  adjustment,  of  the  weighting 
factors,  as  required.  This  keeps  the  number  of  simulated  molecules  more  or 
less  constant  while  allowing  the  number  of  real  molecules  to  adjust  as  the 
solution  evolves.  The  introduction  of  weighting  factors,  with  the  ability 
to  unjust  them  as  the  solution  demands,  is  an  important  feature  ol  a  Monte 
Carlo  simulation  which  is  to  be  usable  by  non-experts. 

4.  GRID  COORDINATES  AND  GRID  STRUCTURE 

As  discussed  in  Section  1,  the  Monte  Carlo  procedure  works  by  breaking 
the  calculation  region  up  into  ceils.  A  solution  cell  should  be  a  region 
in  which  no  properties  change  greatiy,  i.e.,  the  dimensions  of  a  cell 
should  ideally  be  small  compared  to  the  local  scale  length  of  the  flow 
field.  Collisions  are  simulated  on  a  cell-by-cell  basis,  and  molecules  can 
experience  collisions  only  with  other  molecules  in  the  same  cell.  There  is 
no  other  spatial  criterion  used  for  determining  collision  partners,  so  the 
ceil  determines  the  collision  environment  for  any  molecule  within  It.  In 
audition  to  defining  the  collision  environment,  the  other  major  function  of 
the  ceil  structure  is  to  determine  the  points  at  which  output  is  generated. 
There  is  no  requirement  that  the  cells  be  divided  up  in  the  same  coordinate 
system  used  in  the  molecular  state  vector. 

For  Monte  Carlo  calculations,  as  for  other  types  of  computational 
fluid  mechanic  analyses,  the  selection  of  grid  geometry  is  a  critical 
requirement  which  is  often  more  of  an  art  than  a  science.  Considerations 
m  tne  selection  of  a  grid  are: 

•  The  grid  should  be  as  simple  as  possible,  since  the  program 
must  repeatedly  decide  which  ceil  molecules  reside  in  as  they 
move  throughout  the  solution  region.  If  this  determination 
required  the  solution  of  a  complex  equation  or  sifting  through 
tables,  the  entire  program  would  run  significantly  slower  than 
if  the  ceil  can  be  determined  easily. 
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•  The  grid  should  concentrate  cells  where  graaients  are  the 

largest,  so  that  the  least  number  of  total  cells  (and 

molecules)  are  needed  to  obtain  an  accurate  solution. 

•  The  grid  should  provide  flow  field  information  where  it  is 
required,  with  the  resolution  that,  is  desired  for  the  answer 
of  interest. 

'Die  SOCRATES  grid  structure  is  a  simple  extension  of  the  basic 
Cartesian  coordinate  system  that  is  used  elsewhere  in  the  code.  The  cells 
are  determined  by  the  intersection  of  three  families  of  planes,  each  family 
being  perpendicular  to  one  of  the  coordinate  directions.  For  each 

coordinate,  there  is  a  plane  at  zero  and  subsequent  planes  proceed  outward 
in  the  plus  and  minus  coordinate  direction.  For  Instance,  the  intersection 
points,  x j .  on  the  positive  x  axis  are  given  by 


exp( jB/N )  -  1 
XJ  xmax  exp(B)  -  1 


(27) 


where  xm.ju  is  the  position  of  the  last  plane  (the  edge  of  the  solution 
region  ir:  that  direction),  N  is  the  number  of  planes  intersecting  the 
positive  x  axis,  and  B  is  an  adjustable  parameter.  For  B  approaching  zero, 
successive  planes  have  equal  spatial  increments;  and  as  B  is  increased  the 
planes  become  more  concentrated  near  the  origin.  The  same  relation  is 
appiieu  for  planes  intersecting  the  minus  x  direction,  with  xmax  being 
replaced  by  xmjn.  The  other  four  directions  (±y  and  ±z)  are  handled  in  an 
analogous  fashion.  Note  that  the  B  and  N  parameters  are  specified 
separately  for  each  or  the  six  directions  away  from  the  origin,  depending 
on  the  physics  of  the  problem  under  consideration.  These  values  can  be 
input  by  the  user  or  automatically  selected  by  the  program.  The  SOCRATES 
grid  structure  fulfills  the  objectives  enumerated  above  to  a  substantial 
degree.  The  relations  for  cell  boundary  locations  are  easily  inverted  to 
ootuin  the  cell  number  corresponding  to  a  given  position,  and  the  parameter 
ot  the  distribution  allows  for  concentration  of  cells  in  the  inner  region 
wnile  allowing  larger  cells  further  out  where  the  gradients  are  less 
severe . 

A  sample  cell  structure  resulting  from  this  technique  is  illustrated 
in  Figure  3.  For  visual  clarity,  the  number  of  planes  has  been  limited  to 
two  m  each  of  the  six  directions,  since  this  is  the  fewest  number  which 
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Figure  3.  A  Schematic  Showing  the  Basic  Design  of  the  Socrates  Cell 
Structure  Which  Uses  Cartesian  Coordinates  with  Uneven 
Spacing. 

illustrates  the  uneven  spacing.  A  typical  calculation  would  have  several 
times  that  many  planes,  but  the  figure  is  difficult  enough  to  interpret  as 
it  is.  The  shuttle  (or,  in  principle,  any  spacecraft)  can  be  arbitrarily 
located  within  this  cell  structure;  though  it  should  be  located  near  the 
center  at  the  grid  structure  for  it  to  make  sense.  Similarly,  the  wind 
direction  as  seen  from  the  shuttle  can  come  from  any  direction  whatsoever; 
triere  is  nothing  lit  the  ceil  structure  or  coorainate  oefimtion  wnicn 
restricts  it. 


5.  COLLISION  MECHANICS 

5.1  Relations  for  Elastic  Collisions 

The  purpose  of  this  section  is  to  present  relations  appropriate  to  the 
simulation  of  a  collision  in  the  SOCKATES  code.  (The  question  of  how 
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moiocui.es  are  selected  for  collisions,  which  is  cruciai  to  the  proper 
simulation  oi  coilision  frequency,  will  be  taken  up  in  the  next  section.) 
Conservation  of  momentum  implies  that  the  center-of-mass  velocity  of  the 
collision  pair  is  unchanged  by  the  collision;  aria  conservation  of  energy 
then  implies  that  the  magnitude  of  the  relative  velocity  between  the 
collision  partners  is  also  unchanged  ny  the  collision.  Since  the 
collision  is  treated  as  a  statistical  event,  all  that  remains  is  to  select 
tin;  direction  of  the  post -cul lision  relative  velocity  vector  from  the 
correct  distribution.  As  mentioned  in  Section  2,  collisions  in  the  VHS 
mociei  are  treated  as  hard  sphere  collisions  when  they  occur  (though  they  do 
not  occur  with  the  same  velocity  dependence  as  do  hard  sphere  collisions). 
Hence,  as  far  as  the  collision  mechanics  is  concerned,  the  model  is  a  haru 
spnore  model.  For  hard  sphere  molecules,  all  directions  for  the 
post-collision  relative  velocity  vector  are  equally  likely.  Tills  is  the 
chief  computational  simplicity  of  the  VHS  model. 

i.et  the  two  molecules  be  identified  by  subscripts  1  and  2,  with  m  ana 
v  denoting  their  masses  and  velocities.  if  i  and  f  indicate  initial  and 
final  states,  then  the  relations  for  the  collision  can  be  summarized  via: 


*  =  2trg 


(.'12) 


Vincent! ,  W.  G.  and  Kruger,  C.  H  ,  Jr.,  Introduction  to  Physical  Gas 
Dynam i cs ,  John  Wiley  and  Sons,  348  (1965). 
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vr(.cos(e),  sin{e)cos(4>) ,  sin(e)sin(4>)  ] 


(33) 


vt. 


rf 


'if 


Vc,n  +  m. 


m2vrf 


*2 


ana 


'2t 


=  v 


cm 


mlvrf 
n-i  +  m- 


(34) 


(35) 


In  these  relations,  ana  throughout  this  report,  0  indicates  a  ranoom 
variable  which  is  evenly  distributed  on  the  interval  iiero  to  one.  Each 
time  that  0  appears  a  distinct  evaluation  of  the  random  variable  is 
imp. I  leu . 

5.2  Effect  of  Coordinate  System 

Note  that  the  expression  for  the  post-col lision  relative  velocity 
vector  (Eq.  (33))  is  not  coordinate  system  specifie.  The  indicated  vector 
components  can  apply  to  any  locally  orthogonal  coordinate  system,  since  the 
direction  implied  is  random.  The  convenient  coordinate  system  to  use,  of 
course,  is  the  coordinate  system  used  to  define  the  velocity  elements  of 
the  state  vector. 

Although  Eq.  (33)  is  independent,  of  coordinate  system,  there  is  a 
source  of  error  which  is  dependent  on  coordinate  system.  This  error  arises 
from  a  basic  promise  of  the  direct  simulation  Monte  Carlo  methou,  namely 
that  position  in  the  cell  is  ignored  when  selecting  collision  partners.  if 
the  velocities  are  expressed  in  a  coordinate  system  which  has  spatially 
varying  basis  vectors,  then  differences  in  position  between  the  two 
molecules  cart  imply  an  erroneous  difference  in  velocity. 

SOCRATES  makes  use  of  the  effect,  of  spatially  varying  basis  vectors  to 
solve  an  otherwise  difficult  problem.  The  problem  arises  due  to  the 
presence  of  concentrated  sources  of  contaminants,  such  as  thrusters  and 
evaporator  vents,  which  are  modeled  as  point  sources  producing  molecules 
traveling  (initially)  directly  away  from  the  source.  Since  there  is  no 
length  scale  to  a  point  source,  the  assumption  that  properties  are  constant 
Tor  che  ceils  in  the  immediate  vicinity  ol  the  source  must  be  invalid. 
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This  can  result  in  improper  collision  sampling  if  special  care  is  not 
taken . 


It  the  velocities  are  expressed  in  Cartesian  coordinates,  for 
instance,  then  two  molecules  selected  from  different  positions  within  the 
ceil  containing  such  a  point  source  can  have  a  substantial  relative 
velocity.  This  relative  velocity  is  illusory,  however,  since  it  merely 
results  from  the  assumption  of  a  point  source  and  the  neglect  of  spatial 
differences;  there  should  not  be  collisions  based  on  this  relative  velocity 
since  the  molecules  are,  in  fact,  heading  away  from  each  other. 

A  simple  resolution  to  this  problem  is  to  express  the  source  velocity 
elements  in  spherical  polar  coordinates.  In  these  coordinates,  every 
molecule  leaves  the  point  source  with  the  same  velocity  in  the  direction  of 
the  spherical  radius  vector.  Expressed  in  spherical  polar  coordinates,  the 
relative  velocity  disappears.  SOCRATES  transforms  velocity  vectors  to 
spherical  polar  coordinates  for  cells  in  the  vicinity  of  point  sources 
(specifically,  when  the  total  number  density  is  greater  than  three  times 
the  ambient  number  density).  Collisions  are  sampled  in  the  transformed 
coordinates,  and  then  the  velocity  elements  are  transformed  back  to  the 
normal  representation  after  collisions  have  been  sampled  for  the  cell  in 
ques  Hon . 

5.3  Simulation  of  Inelastic  Collisions 

SOCRATES  uses  the  Borgnakke  and  Larsen4  phenomenological  model  for 
transfer  of  energy  between  internal  and  translational  modes.  In  this 
model,  a  collision  is  assumed  to  be  either  perfectly  elastic  or  perfectly 
inelastic,  via  a  user  specified  probability.  A  perfectly  Inelastic 
collision  is  achieved  by  summing  the  total  pre-col lision  energy  (internal 
energy  of  both  molecules  plus  the  translational  energy  of  their  relative 
motion,  Eq.  (4)),  and  then  assigning  post-collision  values  from  the 
equilibrium  distribution  for  collisions  with  that  total  amount  of  energy, 
taking  into  account  the  number  of  internal  degrees  of  freedom  in  the  two 
molecules.  Note  that  this  model  has  the  ability  to  relax  from  a 
nonequilibrium  to  an  equilibrium  state  via  an  effective  collision  number. 
The  ability  to  exchange  internal  energy  in  such  a  manner  comprises  a 
significant  increase  in  capability  for  Monte  Carlo  codes  beyond  the 
previous  models  where  molecules  had  no  internal  energy.  It  is  this 
capability  which  enables  the  codes  to  real  1st ical iy  predict  the  macroscopic 
effects  of  polyatomic  gas  flow. 
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Let  Cj  and  C2  be  the  number  of  internal  degrees  of  freedom  of  the  two 
molecules  in  an  inelastic  collision,  and  Es  be  the  total  collision  energy 
defined  by 

=  £ci  *  £li  +  £2i  1 

where  Eci  is  the  initial  translational  collision  energy  defined  by  Eq.  (4), 
and  Elx  and  £2j  are  the  pre-collision  internal  energies  of  the  two 
molecules.  Using  the  procedures  presented  in  Appendix  A,  the  somewhat 
cumbersome  expressions  given  in  Ref.  (4)  can  be  recast  in  terms  of  the 
chi-square  distribution.  Post-collision  values  for  the  respective  energies 
are  given  by 


Llf 


XiEs 


Xi  +  x2  +  x3 


XoE 


E2f  = 


2*^s 


Xi  +  Xo  +  X, 


(37) 


(38) 


and 


"cf 


x3Es 


Xi  +  X2  +  x3 


(39) 


where  X^  is  selected  from  a  chi-square  distribution  with  Cj  degrees  of 
freedom,  X2  is  selected  from  a  chi-square  distribution  with  C2  degrees  of 
freedom,  ar.o  X3  is  selected  from  a  chi-square  distribution  with  2(2  -  o) 
degrees  of  freedom.  (Efficient  procedures  for  sampling  from  a  chi-squarc 
distribution  are  also  given  in  Appendix  A.)  The  post-coil ision 
translational  energy  is  then  used  to  determine  a  new  relative  velocity 
between  the  two  molecules.  With  this  new  relative  velocity,  the  previous 
relations  for  determining  the  post-collision  velocity  elements  of  the 
molecules  apply  for  inelastic  collisions  as  well  as  for  elastic  collisions. 


The  fact  that  the  translational  energy  Is  selected  from  a  distribution 
with  2(2  -  o)  rather  than  3  degrees  of  freedom  merits  some  explanation.  It 
is  due  to  the  fact  that  these  molecules  are  not  random  samples  from  the 
gas,  but  rather  special  molecules  owing  to  their  being  the  product  of  a 
collision.  This  point  can  perhaps  best  be  seen  by  considering  microscopic 
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reversibility,  where  the  inverse  collision  occurs  with  the  same  rate  in 
equilibrium.  For  this  reverse  process,  molecules  participating  in  it  are 
not  ail  equally  probable,  since  those  with  greater  relative  velocities  are 
more  likely  to  collide.  Hence,  the  number  of  degrees  of  freedom  does  take 
on  tin;  value  three  for  the  special  case  of  o  equal  to  1/2,  which  is 
precisely  the  case  of  collision  frequency  being  independent  of  relative 
velocity.  Translational  energy  in  collisions  behaves  as  if  it  has  2(2  -  o) 
degrees  of  freedom. 

5.4  Collisions  Between  Molecules  with  Distinct  Weighting 
Factors 

There  is  an  obvious  problem  when  considering  a  collision  between  two 
simulated  molecules  with  distinct  weighting  factors,  since  they  represent  a 
different  numoer  of  real  molecules.  If  Wy  and  Wy  represent  the  weighting 
factors  for  the  two  molecules,  with  Wy  being  greater  than  Wy,  then  the 
collision  is  generally  counted  as  Wy  "events".  (More  precisely,  the 
weighting  factor  applied  to  the  collision  is  generally  taken  to  be  Wy.) 
This  is  accomplished  by  always  assigning  post-collision  velocity  and  energy 
components  to  the  state  vector  of  the  molecule  with  the  smaller  weighting 
factor,  but  only  changing  the  components  of  the  molecule  with  the  greater 
we: gating  factor  some  of  the  time.  The  probability  that  the  molecule  with 
tne  greater  weighting  factor  will  have  its  components  changed  is  simply 
Wy/Wy.  Statistically,  this  means  that  for  a  large  number  of  simulated 
collisions,  each  such  simulated  collision  will  average  out  to  Wy  real 
collisions  for  each  species,  even  though  their  weighting  factors  differ, 
it  should  be  noted  that  this  does  violate  conservation  of  momentum  ana 
energy  on  an  individual  collision  basis,  but  these  quantities  are  conserved 
m  the  aggregate  over  a  large  number  of  collisions. 

In  some  cases  the  collision  is  assigned  a  weighting  factor  Wc  which  is 
icss  than  either  of  Wy  or  Wy.  When  this  is  done,  the  velocity  components 
an a  internal  energies  of  the  two  molecules  are  changed  with  a  probability 
of  Wf./Wy  and  Wc/Wy,  respectively.  (See  Section  6.  ) 

5.5  Reactive  Collisions 

Reactive  collisions  can  be  simulated  directly.  The  treatment  of 
reactive  collisions  is  similar  to  that  for  inelastic  collisions,  except 
that  a  heat  of  reaction  is  added  to  the  total  energy  expressed  in  Eq.  (36). 
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Reactive  collisions  can  result  in  the  disappearance  ol  reactant  molecules, 
with  the  post-collision  state  being  applied  to  the  product  molecules. 

5.5.1  Types  of  Reactive  Collisions 

SOCRATES  has  a  fairly  comprehensive  chemistry  package  which  is  capable 
of  handling  a  variety  of  reactive  collisions.  Generally  speaking,  a 
reactive  collision  is  an  event  which  occurs  due  to  collisions  with  a 
probability  that  depends  on  the  velocity  (or  energy)  of  the  collisions. 
The  following  generic  types  of  reactions  are  treatable: 

1.  Specific  Bimolecular  Reactions,  i.e,  reactions  of  the  form 

A  +  B  -»  C  +  D  , 

where  A,  B,  C,  and  D  are  particular  species.  An  example  of 
a  reaction  of  this  type  is 

0  +  H20  -♦  0  +  h2o* 

(In  this  example,  the  vibrationally  excited  state  of  water, 

H20* ,  is  treated  as  a  distinct  species.) 

2.  Generic  Bimolecular  Reactions,  i.e.,  reactions  of  the  form 

A  +  M  -»  B  +  M 

where  A  and  B  are  particular  species,  and  M  can  be  any 
species.  An  example  of  a  reaction  of  this  type  is 

H20  +  M  -»  H2U*  ^  M 

which  is  similar  to  the  previous  reaction  except  that  now 
any  molecule  can  serve  to  excite  the  water  molecule. 

3.  Dissociation  Reactions,  i.e.,  reactions  of  the  form 

A  +  M-*C  +  D  +  M  , 

where  M  is  any  molecule,  and  C  and  D  are  the  fragments  of  A 
that  result  from  dissociation.  An  example  of  a  reaction  of 
this  type  is 

02  +  M-»Q  +  0  +  M 
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5.5.2  Reactive  Coiiision  Probability 


The  Monte  Carlo  program  simulates  all  of  the  above  reaction  types  by 
calculating  a  reactive  cross  section  which  is  a  function  of  the  relative 
collision  energy.  When  a  coiiision  occurs,  the  reaction  is  simulated  with 
a  probability  which  is  proportional  to  the  ratio  of  the  reactive  to 
collision  cross  section  at  the  relative  velocity  for  the  coiiision. 

There  are  two  options  for  specifying  the  reactive  cross  section.  The 
first  is  to  specify  an  Arrhenius  rate  constant,  kr,  of  the  form 

kr  =  ATnexp(-Ea/R0T)  ,  (40) 

where  Ea  is  the  activation  energy  of  the  reaction  and  A  and  n  are 
parameters  of  the  relation.  (Rq  is  the  gas  constant  and  T  is  temperature.) 
The  unique  reactive  cross  section,  <r  ,  corresponding  to  Eq.  (40)  is  given 
by3 


(1  t  SijJtr-SA 
2R0nP(n  +  3/2) 


vra 


t/1  -  Ep./Ef,  (hn  -  Ea)n 


(41) 


where  Sjj  is  unity  for  like  reactants  and  zero  for  unlike  reactants,  T 
represents  the  gamma  function  and  Ec  is  the  collision  energy  given  by  Eq. 
(4).  Note  that  a  rate  constant  is  defined  in  terms  of  an  equilibrium 
velocity  distribution,  so  the  correspondence  between  Eqs.  (40)  and  (41)  can 
be  maoe.  There  is  no  requirement,  of  course,  that  the  reactive  cross 
section  given  by  Eq.  (41)  be  used  only  in  equilibrium  situations.  When 
this  option  is  used,  only  the  arrhenius  parameters  A,n  and  Ea  need  be 
specified;  the  program  automatically  computes  the  corresponding  reactive 
cross  section. 


Eor  some  reactions,  the  form  of  Eq.  (41)  is  too  restrictive,  and  it  is 
then  possible  to  input  a  table  giving  the  reaction  cross  section.  The  form 
of  the  table  is  of  the  same  functional  form  as  Eq.  (41),  namely  the  product 
of  the  relative  velocity  times  the  reactive  cross  section  is  given  as  a 
function  of  relative  collision  energy.  Although  this  form  is  not  standard, 
■t  .s  far  more  convenient  for  reactions  where  one  of  the  reactants  is 
generic  ( "M " ) ,  since  there  is  no  correspondence  between  collision  velocity 
and  collision  energy  until  the  masses  of  both  reactants  are  specified. 
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5.5.3  Options  for  Simulating  Reactive  Collisions 

SOCRATES  has  distinct  options  for  simulating  reactive  collisions  which 
are  reflective  of  different  anticipated  user  needs.  In  all  options,  the 
sampling  of  the  reaction  rate  (if  it  is  being  performed)  is  done  the  same 
way.  Whenever  a  collision  occurs  between  the  two  reactants,  the  reactive 
cross  section  is  calculated,  and  the  reaction  is  counted  with  a  weighting 
factor,  Wr,  given  by 

* 

vrff 

Wr  =  Wc  — .  (42) 

r  c  vrff 

where  Wc  is  the  collision  weighting  factor  (see  the  previous  section). 
Hence,  even  though  the  reactive  cross  section  may  be  signif icantly  smaller 
than  the  collision  cross  section,  the  statistics  on  the  reaction  rate  are 
similar.  (The  statistics  for  the  reaction  rate  may  converge  slower  due  to 
the  velocity  dependence  of  the  reaction  cross  section;  but  not  due  to  its 
absolute  magnitude.)  If  two  molecules  can  participate  in  multiple 
reactions,  statistics  are  kept  for  each  reaction. 

If  products  are  introduced  as  a  result  of  the  reaction,  they  can  be 
introduced  at  every  simulated  reactive  collision  with  a  weighting  factor  of 
Wr,  or  introduced  with  a  weighting  factor  of  Wc,  but  only  Wr/Wc  of  the 
time.  The  difference  depends  on  the  importance  of  tracing  product  species 
in  the  simulation.  The  former  approach  will  result  in  more  computational 
effort  being  spent  on  the  product  species,  but  it  will  give  better 
statistics  on  them.  In  either  case,  reactants  are  removed  from  the 
simulation  with  a  probability  of  Wr/Wc  in  any  reactive  collision. 

In  many  cases,  it  is  the  reaction  rate  which  is  of  interest.  If  the 
reactive  collisions  are  relatively  improbable  events,  then  the  velocity 
distribution  of  the  reactants  will  be  the  same  irrespective  of  whether  the 
reaction  is  actually  simulated.  For  T-V  excitation  reactions,  which  are  of 
primary  interest  for  IR  interference,  it  is  possible  to  calculate  the 
excitation  rate  without  explicitly  introducing  excited  state  molecules  into 
the  simulation,  or  removing  the  ground  state  molecule  which  becomes 
excited.  The  calculation  of  reaction  rates  in  such  a  manner  is  easily 
incorporated  into  both  procedures  for  simulating  equilibrium  flow. 


Furthermore,  such  a  procedure  can  be  extended  to  give  accurate 
estimates  for  the  competing  processes  of  collisional  quenching  and 
radiative  decay.  Since  a  vi brationally  excited  molecule  has  the  same  mass 
as  its  ground  state  counterpart,  the  velocity  distribution  of  excited  state 
molecules  should  be  well  approximated  by  velocity  distribution  of  the 
corresponding  ground  state  molecules.  Hence,  it  is  possible  to  define  an 
artificial  quenching  reaction  which  has  ground  state  molecules  as  both 
reactants  and  products  (but  with  the  proper  quenching  cross  section).  To 
make  the  discussion  more  concrete,  consider  the  following  reaction  scheme: 


A  +  M  -+  A  +  M 


A  +  M  -*  A  +  M 
,  kr 

A  -♦  A  +  hv 


(Rl) 


(R2) 


(R3) 


where  excitation,  quenching,  and  radiative  decay  are  represented  by 
Reactions  1  through  3,  respectively.  The  rate  constants  ke  and  kq  are 
determined  by  the  existing  velocity  distribution  of  the  reactants  and  the 
relevant  cross  sections  as  a  function  of  collision  velocity.  (The  rate 
constant  kr  is  simply  the  inverse  of  the  radiative  lifetime.  Tr.)  Hence, 
if  the  velocity  distribution  of  the  excited  state  species  is  the  same  as 
for  the  ground  state,  then  the  rate  constants  are  the  same  in  the 
artificial  reaction  scheme: 


At+  M 


*'e 

-*  A. 


+  M 


(Rl ') 


A^  +  M 


At+  M 


(R2*  ) 


wnere  At  now  represents  both  ground  state  and  excited  state  molecules.  The 
artificial  reaction  scheme  does  not  require  the  removal  of  reactants  or  the 
introduction  of  products,  but  it  does  provide  the  basic  rate  constant 
information  (i.e.,  ke  and  kq)  required  to  determine  emission  rates.  At 
equilibrium,  Reactions  2  and  3  will  balance  Reaction  1,  and  A  will  be 
re lat  :d  to  A  by 
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(43) 


Ake  =  A*(kq  +  kr) 

Substituting  the  definition  of  At , 

At  =  A  +  A*  (44) 

it  is  possible  to  solve  for  the  emission  rate  via 


A*kr 


kekr 

ke+kq+krAt 


(45) 


Ail  of  the  information  required  to  calculate  the  emission  rate  via  Eq.  (45) 
is  obtainable  via  a  chemical  reaction  scheme  in  which  products  are  never 
explicitly  introduced  and  reactants  are  never  explicitly  removed.  This 
capability  is  an  option  in  SOCRATES,  and  it  is  completely  compatible  with 
the  procedures  for  speeding  up  execution  in  collision  dominated  scenarios. 


6.  COLLISION  SAMPLING  IN  A  MULTI-COMPONENT  VHS  GAS 
6.1  General  Considerations  and  Approach 

The  two  general  considerations  in  the  sampling  of  collisions  are,  as 
usual,  accuracy  and  efficiency  of  the  simulation.  As  far  as  accuracy  is 
concerned,  it  is  crucial  that  the  method  in  which  molecules  are  selected 
for  collisions  be  proper.  The  correct  collision  frequency  must  be 
simulated  between  various  species  and,  in  fact,  between  the  different 
portions  of  the  velocity  phase  space  for  the  various  species.  Furthermore, 
this  frequency  of  simulated  collisions  must  remain  correct  without  any 
requirements  put  on  the  velocity  distribution  function;  it  certainly  must 
not  be  assumed  that  there  is  a  Maxwellian  velocity  distribution. 

As  far  as  efficiency  is  concerned,  it  is  highly  desirable  to  use  a 
method  of  collision  sampling  involving  a  computational  effort  which  is 
proportional  to  the  number  of  simulated  molecules,  N,  in  a  cell.  Methods 
which  are  proportional  to  a  power  of  N  greater  than  unity  can  become 
prohibitively  time  consuming  as  the  number  of  molecules  is  increased  -  a 
limit  which  should  be  made  as  accessible  as  possible  for  obvious  physical 
reasons . 
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6.2  Collision  Sampling  for  a  Single  Component  Gas 

the  simplified  situation  of  a  simulation  involving  only  one  species  is 
considered  here.  This  problem  is  significant  in  part  due  to  all  the 
attention  it  has  received  and,  as  will  be  seen,  it  serves  as  an  important 
reference  case.  When  there  is  just  one  species,  then  there  is  just  one  gas 
kinetic  cross  section  (though  it  is  still,  of  course,  a  function  of 
collision  energy),  just  one  molecular  weight  and  just  one  weighting  factor 
for  each  cell.  In  short,  just  one  of  everything  that  has  a  molecular 
subscript.  Hence,  in  this  subsection  all  such  quantities  will  be  presented 
without  subscripts.  The  most  important  simplif ication  of  having  a  single 
species  is  that  there  is  just  one  collision  class,  i.e.,  only 
sui f -col  1 isions  of  the  given  species  with  itself  are  possible. 

b.2.1  Collision  Pair  Selection 

As  discussed  in  Section  1,  collisions  are  sampled  on  a  cell-by-cell 
basis  until  the  number  of  collisions  simulated  is  appropriate  to  the 
overall  solution  time  step,  htm.  The  only  spatial  requirement  placed  on 
potential  collision  partners  is  that  they  be  within  the  same  cell.  In 
particular,  it  is  not  required  that  they  be  within  a  molecular  diameter  of 
each  other.  (Note  that  if  all  pairs  of  molecules  were  inspected  to  find 
those  that  were  sufficiently  close  to  each  other,  this  would  involve  a 
computational  effort  in  proportion  to  the  square  of  the  number  of  molecules 
in  the  cell.)  The  rationale  for  this  is  that  the  cells  should  be  small 
enough  so  that  macroscopic  properties  can  be  assumed  constant  across  the 
ceil.  When  this  is  the  case,  then  a  molecule  within  the  cell  can  be 
considered  typical  of  a  molecule  which  might  exist  anywhere  within  the 
cell,  and  molecular  location  can  be  ignored  when  selecting  potential 
collision  pairs. 

Spatial  consideration  aside,  the  probability  of  any  two  molecules 
experiencing  a  collision  is  proportional  to  <rcr,  the  product  of  their 
mutual  cross  section  times  their  relative  velocity.  This  probability  is 
correctly  simulated  via  an  application  of  the  acceptance-rejection 
technique.  A  maximum  value  for  trcr,  is  stored  for  each  cell. 
(This  value  is  updated  whenever  a  greater  value  is  encountered.)  Pairs  of 
molecules  are  selected  at  random  from  the  cell,  and  <rcr  for  that  pair  is 
calculated.  The  ratio  r,  defined  by 
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(ffcr  ^max 

is  determined.  A  random  variable,  0,  is  then  generated,  and  the  pair  of 
molecules  is  accepted  as  collision  partners  if  r  is  greater  than  0.  This 
produces  the  proper  relative  collision  probability  without  regard  to  the 
existing  velocity  distribution  function. 

6.2.2  Collision  Time  Counter  for  a  Single  Component  Gas 

The  volumetric  collision  frequency  for  a  single  component  gas,  v, 
(collisions  per  unit  volume  per  unit  time),  is  given  by 
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where,  as  in  Section  2,  n  represents  the  number  density  of  the  species,  and 
<ffcr>  is  the  average  product  of  collision  cross  section  and  relative 
velocity.  At  first  inspection,  it  would  seem  from  Eq.  (47)  that  a  correct 
simulation  of  collision  frequency  would  require  evaluation  of  <ffcr>,  which 
would  mean  that  all  pairs  of  molecules  in  a  cell  would  have  to  be 
considered.  Such  a  procedure  involves  a  computational  effort  proportional 
to  and  is  to  be  avoided,  if  possible,  in  preference  to  a  method  which  is 
proportional  simply  to  N.  The  alternative  approach,  introduced  by  Bird,1 
is  the  time  counter  approach.  For  each  collision  a  time  counter,  tc ,  is 
incremented  by  an  amount  which  depends  on  the  relative  velocity  of  the 
collision.  Collision  sampling  continues  in  a  cell  until  its  time  counter 
has  been  advanced  beyond  the  overall  flow  simulation  time,  at  which  time 
the  code  proceeds  to  the  next  cell  (which  has  its  own  time  counter).  The 
time  counter  increment,  Atc,  is  given  by 
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where  V  is  the  cell  volume  and  n  is  the  species  number  density  given  by 
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with  W  being  the  weighting  factor  for  the  species.  It  should  be  stressed 
that  Eq.  (48)  applies  for  each  real  collision.  As  is  discussed  in 
Subsections  3.4  and  5.4.  each  simulated  collision  corresponds  to  W  real 
collisions,  so  when  a  simulated  collision  occurs  the  actual  applied 
increment,  to  cc  is  W  times  the  value  given  by  Eq.  (48).  A  demonstration  ol 
the  validity  of  Eq.  (48)  is  given  in  Ref.  6. 

6.3  Collision  Class  Sampling  in  Gas  Mixtures 

The  above  procedure  for  a  single  species  gas  can  be  extended  to  a 
multi-component  mixture  via  consideration  of  distinct  collision  classes. 
In  this  approach,  collision  classes  are  defined  by  the  colliding  pair 
identities.  Hence,  if  there  are  p  species  in  the  simulation  then  there  are 
p(p-»l)/2  collision  classes,  which  can  be  identified  by  the  subscripts  of 
the  corresponding  molecular  pair.  (The  number  of  classes  is  not  p2  since 
the  oroer  of  molecule  specification  is  not  taken  to  matter  in  determining  a 
collision  class.  Hence,  the  class  identified  by  the  subscripts  i,j  is  not 
distinct  from  the  class  identified  by  the  subscripts  j.i.) 


In  collision  class  sampling  each  collision  class  is  sampled 
separately,  and  the  collision  sampling  in  a  cell  is  not  complete  until  all 
classes  have  been  considered.  Each  collision  class  has  its  own  stored 
value  of  (o'ijCr)max  and  its  own  separate  time  counter,  tci  j .  It  can  be 
shown  that  the  appropriate  time  counter  increment  in  this  case  is 
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wnere,  as  before,  S^j  is  the  Kronecker  delta  which  is  unity  for  i*j  and 
zero  otherwise.  As  in  the  previous  section,  the  above  increment  applies 
tor  each  real  collision.  A  simulated  collision  usually  corresponds  to  WL 
real  collisions,  where  is  the  lesser  of  and  Wj  (see  Subsection  5.4), 
so  when  a  simulated  collision  occurs,  the  applied  increment  to  tc^j  is 
times  the  result  of  Eq .  (50). 


6  Elgin,  J.  B. ,  "Getting  the  Good  Bounce:  Techniques  for  Efficient 

Monte  Carlo  Analysis  of  Complex  Reacting  Flows',  Spectral  Sciences, 
Inc.  Report.  No.  SSI-TK-28  (1983). 


29 


6.4  Global  Collision  Sampling  in  a  Gas  Mixture 


AJ though  the  procedure  described  above  is  quite  reasonable  for,  say,  a 
two-component  mixture,  it  becomes  exceedingly  complicated  as  the  number  of 
species  increases.  For  10  species,  for  instance,  the  program  must  loop 
over  55  distinct  collision  classes  for  each  cell,  and  storage  must  be 
ai located  for  110  quantities  in  each  cell.  As  the  number  of  species 
increases,  the  storage  requirement  for  the  collision  sampling  constants 
quickly  becomes  greater  than  the  storage  required  for  the  molecular  state 
vectors!  The  obvious  simplification  is  to  search  for  a  technique  where 
collisions  are  simuiated  simultaneously  for  all  collision  classes,  with 
each  class  having  its  proper  relative  probability  of  being  selected.  The 
overall  collision  sampling  then  continues  until  a  single  time  counter 
indicates  that  sufficient  collisions  have  been  sampled  in  the  current  time 
step  ana  cell. 

6.4.1  Global  Collision  Time  Counter 

If  molecular  pairs  are  selected  for  collisions  such  that  the  various 
collision  classes  automatically  appear  with  the  proper  relative  frequency 
(see  below),  then  it  is  not  necessary  to  consider  separate  time  counters 
for  ail  the  various  collision  classes.  One  approach  that  could  then  be 
applied  is  to  i  eep  a  collision  time  counter  for  just  one  collision  class 
and  increment  it  when  collisions  of  that  class  occur.  If  the  various 
collision  classes  are  being  selected  according  to  their  correct  relative 
frequency,  tnen  simulating  the  proper  frequency  for  one  collision  class 
will  ensure,  in  the  long  run,  that  ail  collision  classes  are  occurring  with 
the  correct  frequency.  A  disadvantage  with  this  approach  is  the  necessity 
of  making  an  arbitrary  choice  for  the  collision  class  which  is  to  have  a 
time  counter.  Furthermore,  there  may  be  no  good  choice  for  a  reacting  flow 
where  the  dominant  species  can  vary  strongly  from  place  to  place. 
(Clearly,  one  would  not.  want,  to  select  a  class  of  collision  that,  does  not 
occur  in  u  given  cell,  since  the  result  would  be  a  never-ending  samp'ing  of 
collisions  of  other  classes.) 

The  preferred  approach  is  to  define  a  global  collision  time  counter, 
t  which  is  a  weighted  average  of  the  time  counters  of  all  collision 
classes ;  i . e . , 
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(51  ) 


where 


P  i 

c  =  y  yD^  ,  (52) 

,]=i  j=i 

and  the  Djj  are  non-negative  coefficients  which  can  be  selected  at  will. 
Nuie  that  in  this  formulation  every  collision  class  will  result  in  some 
increment  of  the  global  time  counter  (unless  D j j  is  zero  for  that  class), 
so  the  collision  sampling  frequency  is  not  dependent  on  any  one  collision 
class . 

it  remains,  of  course,  to  specify  the  D^j.  A  very  convenient  choice 
is  given  by 


(53) 


Firstly,  Eq.  (53)  is  convenient  because  it  tends  to  make  the  collision 
classes  with  the  higher  collision  frequencies  count  more,  resulting  in  goon 
statistics  for  tg  irrespective  of  cell  location.  (Note  that  is  cell 
dependent  since  the  species  number  densities  are  ceil  dependent.) 
Secondly,  Eq.  (53)  results  in  a  particularly  convenient  form  for  tg.  The 
normalization  factor  given  in  Eq.  (52)  can  be  summed  analytically  to  give 


P  i 

■c-'r  v*  nini 

>Z 

i=l  j=l 


Hence,  a  collision  of  class  ij,  which  wouid  produce  an  increment  of  Atcij 
to  its  own  time  counter  produces  an  increment  Atg  to  tg  given  by 


where,  again,  n  is  the  total  number  density  of  all  species  in  the  cell.  If 
Eq.  (50)  is  substituted  into  Eq.  (55),  the  result  is 
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Equation  (56)  is  extremely  significant  since  it  recaptures  the  precise  form 
of  the  time  counter  increment  for  a  single  species  (Eq.  (48)),  but 
indicates  that  it  is  completely  valid  for  a  mul ti -component  mixture  so  long 
as  the  various  collision  classes  are  sampled  with  the  proper  relative 
frequency . 


6.4.2  Collision  Fair  Selection  in  Multi-Component  Mixtures 


When  considering  selection  of  collision  pairs,  it  is  crucial  to 
remember  the  distinction  between  real  and  simulated  molecules  discussed  in 
Subsection  3.4.  Given  two  simulated  molecules  selected  at  random  from 
within  the  ceil,  the  probability  of  their  having  a  real  collision  is 
proportional  to  W|WjffjjCr.  However,  real  collisions  cannot  happen 
individually;  they  come  WL  at  a  time,  where  WL  is  the  lesser  of  and  W^. 
Hence,  wnen  a  collision  is  decided  upon  in  the  program,  WL  of  them  will 
occur.  To  compensate  for  this,  potential  collision  pairs  should  be 
accepted  for  a  collision  according  to  the  size  of  Q  given  by 


Q  15  wUffij(:r  (57) 

The  reiative  frequency  of  real  ij  collisions  will  then  be  proportional  to 
the  product  QW^  (the  relative  probability  of  a  pair  being  accepted  for  a 
collision  times  the  number  of  real  collisions  occurring  when  the  pair  is 
accepted),  which  is  the  desired  relation.  Selection  of  collision  pairs 
with  the  correct  relative  frequency  then  assures  that  incrementing  the 
global  time  counter  as  discussed  above  will  give  a  statistical ly  correct 
sampling  of  ail  collision  classes  simultaneously. 

6.4.3  Summary  of  Collision  Sampling  in  Multi-Component  Mixtures 

The  results  of  this  subsection  can  be  summarized  via  the  following 
procedure  for  the  sampling  of  collisions: 

•  Each  cell  has  a  (current)  maximum  value  of  Q,  Qmax.  that  has 
been  encountered  so  far  in  the  collision  sampling  process. 
Whenever  a  larger  value  is  encountered,  Qmax  is  set  equal  to 
that  larger  value. 
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•  Each  cell  has  a  current  value  of  the  global  time  counter,  tg. 

•  Pairs  of  simulated  molecules  are  selected  at  random  from  ail 
molecules  within  the  cell. 

•  For  each  pair,  Q,  (as  defined  by  Eq.  (57))  is  computed. 

•  The  ratio  of  Q  to  Qmax  is  computed,  and  a  random  variable  is 

generated.  The  pair  is  accepted  for  collision  if  the  random 
variable  is  less  than  that  ratio.  (If  the  pair  is  not 

acceptea,  then  anocher  random  pair  is  selected.  The  process 
continues  until  a  pair  is  accepted.  ) 

•  For  an  accepted  pair,  the  collision  mechanics  are  computed  as 
described  in  Section  5. 

•  The  global  time  counter  is  incremented  by  W^dtg,  where  Atg  is 
given  in  F,q.  (56). 

•  The  process  continues  until  the  global  time  counter  goes 

oeyond  the  overall  flow  time.  At  that  point  the  collision 
sampling  is  commenced  in  the  next  cell. 

•  When  all  cells  have  had  collisions  simulated,  then  the  code 

proceeds  to  the  translation  portion.  (See  Sections  1  and  8.) 

6.5  Deviations  from  the  General  Procedure 

There  are  some  exceptions  to  the  above  relations  which  have  been  added 
to  SOOKATES  in  order  make  it  more  efficient.  These  exceptions  are 

described  in  the  following  subsections. 

6.5.1  Cell  Specific  btm 

Before  collisions  are  simulated  in  a  cell,  the  mean  residence  time  of 
molecules  in  the  cell  is  estimated  using  the  cell  dimensions  and  the 
molecular  velocities.  When  collisions  are  simulated  in  the  cell,  it  is 

done  for  an  increment  of  the  global  time  counter  that  is  20%  of  this  mean 

residence  time  (but  no  less  than  Atm).  Collisions  are  not  again  simulated 
in  the  ceil  until  the  overall  flow  time  has  caught  up  to  the  global  time 
counter  tor  the  cell. 

The  major  reason  for  doing  this  is  to  recognize  that  some  cells  will 
tend  to  have  their  molecules  remain  in  them  much  longer  than  others.  Cells 
which  have  longer  molecular  residence  times  will  tend  to  have  molecules 
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which  experience  more  collisions  within  the  cell.  When  the  number  of 
collisions  per  molecule  becomes  sufficiently  large,  it  can  be  assumed  that 
the  molecules  in  the  cell  equilibrate  with  each  other,  and  the  equilibrium 
sampling  procedures  described  in  Section  7  can  be  applied.  Since  hese 
relations  are  much  faster  than  direct  collision  sampling,  it  is  highly 
desirable  to  apply  them  whenever  they  are  valid. 

For  unsteady  simulations  the  cell  specific  Atm  is  not  applied  since  it 
mignt  result  in  a  temporal  blurring  of  the  solution. 

6.5.2  Relaxation  of  Qmax 

The  current  value  of  Qmax  in  a  cell  is  reduced  by  a  factor  of  0.95  if 
20  or  more  potential  collision  pairs  are  rejected  in  a  row.  The  rejection 
of  collision  pairs  can  become  the  most  time-consuming  part  of  the 
simulation,  and  a  large  value  of  Qmax  exacerbates  the  problem.  This  change 
means  that  a  cell  is  not  permanently  penalized  for  a  single  event  that  once 
occurred  in  it,  but  the  change  in  Qmax  is  not  so  great  as  to  invalidate  the 
pair  selection  probability.  This  modification  can,  under  some  collision- 
aominated  circumstances,  result  in  an  order  of  magnitude  increase  in 
computational  speed. 

6.5.3  Maximum  Time  Counter  Increment 

Since  Atg  is  inversely  proportional  to  relative  velocity  (Eq.  (56)), 
when  a  very  low  velocity  collision  does  occur,  it  can  result  in  very  large 
increment  to  the  collision  time  counter,  which  effectively  turns  off 
collisions  in  the  ceil  for  a  long  time.  Although  this  Is  statistically 
proper  in  the  long  run.  It  can  result  in  a  substantial  statistical 
fluctuation  in  the  short  run.  The  codes  do  not  allow  a  collision  time 
increment  to  be  greater  that  htm,  the  overall  step  that  is  used  in  the 
solution . 

The  limitation  on  Atg  is  achieved  by  decreasing  the  weighting  factor 
ol'  the  collision  below  the  weighting  factor  of  either  of  the  two  colliding 
molecules.  The  maximum  collision  weighting  factor,  (Wc)max.  is  given  by 
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The  weighting  factor  that  is  applied  to  a  collision  is,  therefore,  the 
smaller  of  <Wit  Wj  ,  (wc)max*-  wc  represents  this  value,  then  a 

collision  counts  as  Wc  events.  In  order  to  maintain  an  overall  correct 
simulation,  the  Q  described  in  Eq .  (57)  is  actually  given  by  the  relation 
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ana  the  time  increment  applied  to  the  global  time  counter  is  W(.Atg.  (The 
two  molecules  then  have  their  state  vectors  updated  as  a  result  of  the 
collision  with  probabilities  of  Vlc/V^  and  Wc/Wj  respectively.)  Most  of  the 
time  Wc  is  equal  to  W^,  and  this  procedure  reduces  to  that  given  above; 
however,  the  problem  of  occasional  large  time  increments  is  eliminated. 

6.5.4  Separation  of  Major  and  Minor  Species 

A  problem  arose  when  a  cell  happened  to  contain  a  single  molecule  of 
one  major  species  and  all  other  molecules  in  the  cell  were  minor  species 
with  weighting  factors  considerably  less  than  that  of  the  first  molecule. 
(In  treating  minor  species,  the  ratio  of  weighting  factors  may  be  as  much 
as  1000  or  more.)  Since  a  moiecuie  cannot  collide  with  itself,  collisions 
between  major  species  could  not  occur  in  such  a  cell.  The  result  was  that 
the  contribution  of  major  species  collisions  to  the  overall  time  counter 
were  unobtainable;  ana  the  entire  collision  time  increment  had  to  be  made 
up  with  collisions  between  the  single  major  species  and  one  or  other  of  the 
minor  species.  (Collisions  between  minor  species  were  rare  since  pairs  are 
selected  with  a  probability  which  is  proportional  to  the  greater  weighting 
factor  -  see  above.)  The  result  of  this  problem  was  that  vastly  too  many 
collisions  were  simulated  between  the  major  and  minor  species.  This  was 
both  unpnysicai  and  numerically  inefficient. 

The  solution  to  the  problem  is  twofold:  1)  Logic  is  in  the  modules  to 
recognize  when  this  problem  occurs;  and  2)  The  global  time  counter  is 
redefined  in  such  situations.  The  different  time  counter  is  achieved 
simply  choosing  a  different  definition  for  the  Djj  coefficients  appearing 
in  Eq.  (51).  Rather  than  taking  a  weighted  average  over  all  collision 
classes,  it  is  possible  to  take  a  weighted  average  over  just  those 
collision  classes  which  involve  a  collision  between  a  major  and  a  minor 
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species.  It  n^  represents  the  total  major  species  number  density  and  nj 
represents  the  total  minor  species  number  density,  then  Djj  is  defined  for 
this  case  to  be 

Dij  =  ninj  -  (60) 

rather  than  the  value  given  in  Eq.  (53).  The  implied  increment  for  the 
global  time  counter  is  given  by 
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This  counter  is  then  only  applied  for  collisions  between  the  species 
declared  to  be  major  and  the  species  declared  to  be  minor,  but  the 
increment  that  is  applied  is  much  larger  than  if  the  weighting  were  over 
all  collision  classes.  Note  that  collisions  between  minor  species  still 
occur  -  they  just  don't  affect  the  collision  time  counter.  It  should  be 
stressed  that  this  modification  is  only  applied  for  the  special  case  of  a 
cell  that  has  a  single  major  species  molecule  (defined  as  a  weighting 
factor  at  least  ten  times  greater  than  that  of  the  other  molecules).  The 
result  of  this  modification  was  the  return  of  the  correct  collision 
frequency  simulation  and  a  removal  of  a  substantial  numerical  problem. 


7.  PROCEDURES  FOR  COLLISION  DOMINATED  FLOW 

One  of  the  major  difficulties  in  the  classical  Monte  Carlo  technique 
is  the  attainment  of  equilibrium,  where  the  collision  frequency  can  become 
prohibitively  large  for  a  direct  simulation.  There  are  two  basic 
approaches  tor  this  problem,  and  both  are  utilized  in  SOCRATES.  The 
equilibrium  modeling  is  only  applied  if  no  products  are  being  introduced 
info  the  simulation  as  a  result  of  chemical  reactions.  When  such  prouucts 
are  introduced,  collisions  are  always  sampled  for  the  full  time  increment. 
The  two  equilibrium  approaches  are  described  below. 

7.1  Collision  Cutoff  Approach 

This  is  the  usual  method  of  dealing  with  a  collision  dominated  flow 
field.  In  this  method,  collisions  are  sampled  in  a  given  celi  only  until 
enough  have  been  sampled  to  guarantee  equilibrium.  Since  further 
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collisions  only  result  in  the  maintaining  of  equilibrium,  they  need  not  be 
simulated.  It  is  necessary,  of  course,  to  estimate  the  actual  collision 
frequency  and  keep  proper  count  of  the  collisions  (in  particular  the 
excitations)  which  are  not  directly  simulated  in  order  to  obtain  the 
correct,  collision  frequency.  Once  a  ceil  has  had  its  collisions  cut  off  in 
this  fashion,  a  flag  is  set  for  it.  On  subsequent  calls,  the  equilibrium 
aftermath  approach  (described  in  the  next  section)  is  applied  to  the  cell. 
(The  equilibrium  aftermath  approach  also  calculates  collision  frequencies 
ana  switches  back  to  regular  collision  sampling  if  it  becomes  too  low.) 

7.2  Equilibrium  Aftermath  Approach 

It  is  possible  to  avoid  sampling  collisions  altogether  if  it  is  known 
that  the  cell  is  in  equilibrium.  This  is  done  by  calculating  the  total 
cell  energy  and  momentum,  and  then  selecting  post-collision  velocities  from 
the  appropriate  equilibrium  distribution.  Although  the  principle  is 
simple,  the  application  is  complicated  by  the  fact  that  molecules  in  the 
ceil  do  not  all  have  the  same  statistical  weight.  In  some  ways  (e.g.,  the 
determination  of  mean  velocity)  the  statistical  weight  acts  like  an 
effective  multiplication  of  molecular  mass  -  the  greater  a  molecule's 
statistical  weight,  the  greater  its  contribution  to  the  mean  flow  velocity. 
In  other  ways  (e.g.,  the  assignment  at  post-collision  thermal  velocities) 
the  statistical  weight  does  not  affect  the  result  -  a  light  molecule  should 
generally  have  a  large  thermal  velocity  irrespective  of  its  statistical 
weight . 


The  second  difficulty  with  the  formulation  of  this  approach  is  the 
necessity  for  constraining  the  total  energy  and  momentum  to  match  the 
pre-collision  values.  Hence,  it  is  not  proper  to  calculate  the  initial 
energy  and  then  simply  sample  from  Maxwellian  distributions  with  the  same 
mean  energy,  since  such  a  distribution  has  a  finite  probability  of 
producing  a  molecule  with  any  energy  -  and  the  net  result  would  be  an 
unacceptable  divergence  in  the  cell  energy  and  momentum  from  the  initial 
values.  Both  of  these  problems  are  avoided  in  the  steps  enumerated  below. 

The  method  is  implemented  by  calculating  the  total  momentum  and  energy 
in  the  ceil  and  then  "peeling  off"  one  molecule  at  a  time  from  the  others. 
The  .infernal  mode  energy  and  energy  of  relative  motion  for  that  molecule 
(relative  motion  with  respect  to  the  remaining  molecules)  are  selected  from 
equilibrium  distributions,  except  that  a  scale  factor  which  is  proportional 
to  temperature  is  temporarily  left  undetermined.  The  process  is  repeated 
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seouentially  until  all  energy  mones  have  been  assigned  values,  and  then  the 
overall  multiplicative  constant  is  chosen  to  match  the  known  total  energy 
ol  the  system  (thus  determining  the  temperature). 

Each  molecule  is  then  assigned  velocity  components  which 
consistent  with  the  known  relative  velocity  between  it  and  the  remaining 
molecules;  and  then  conservation  of  momentum  determines  the  mean  velocity 
of  the  remaining  molecules.  Again,  the  process  is  repeated  sequentially 
until  all  velocities  and  internal  energies  have  been  assigned. 

7.2.1  Conserved  Quantities 

The  total  energy  and  center-of-mass  velocity  are  directly  computed  via 
the  following  procedure: 

1)  The  following  sums  are  evaluated,  summing  over  ail  the  simulated 

molecules  in  the  cell: 


S1  “  ^wimi 

s2  =  Z«i»iui  • 

S3  •  2"i»ivi  • 

s4  =  2wimiwi 

s5  -  7  Wjfflj  (u“-  +  v*  +■  w^> 
arid 

s6  -  - 


(62) 

(63) 

(64) 

(65) 

(66) 

(67) 
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where  Wj ,  nij  and  Ejj  are  the  statistical  weighting  factor,  the 
mass  ana  the  internal  energy,  respectively,  of  the  ith  molecule; 
ana  U| ,  v^  and  wa  are  its  velocity  components. 

$  £  ♦ 

2)  The  center  of  mass  velocity  components,  u  .  v  and  w  ,  are 
computed  via: 

u*  =  Sg/Sj  ,  (68) 

v"*  =  S3/S1  (69) 

and 

w*  =  84/8!  (70) 

3)  The  total  translational  energy  of  relative  motion  between  the 
molecules,  Etrn,  can  be  represented  by: 


=  |ywiiai[(ui  -  u*)2  +  (v^  -  v*)2  +  (w^  -  w*)2] 


although  it  is  more  easily  evaluated  by  the  mathematically 
equivalent  expression 


E 


trn 


(72) 


4)  The  total  cell  energy  is  therefore  given  by 


Etot  =  s6  +  Etrn 


(73) 


7.2.2  Center-of-Mass  Velocity  Distribution 

Given  that  a  group  of  N  molecules  is  in  equilibrium,  it  is  possible  to 
determine  the  form  of  the  distribution  function  for  their  mass  averaged 
velocity,  taking  into  account,  their  different  statistical  weighting 
factors.  This  relation  is  most  easily  demonstrated  by  relating  the 
Maxwellian  velocity  distribution  to  the  normal  distribution  of  statistics 
ana  then  utilizing  a  basic  statistical  theorem. 
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A  variable,  r,  is  distributed  according  to  a  normal  distribution  it 
its  probability  density  function,  f(r),  is  given  by 


t*2  _ 

f(r)  =  exp(-~-)y  {a\'2rr )  ,  (74) 

2<r2 

where  <r2  is  the  variance  of  the  distribution.  (The  distribution  has  been 

selected  with  zero  mean  since  the  effect  of  non-zero  means  does  not 

influence  the  velocity  differences  which  are  the  goal  of  this  exercise.)  A 

basic  result  of  statistics  is  that  if  r,  is  selected  from  a  normal 

2  x 

uistribution  with  variance  O',  and  r->  is  selected  from  a  normal  distribution 
2  16 

with  variance  a,,,  then  the  variable  r3  defined  by 

r3  =  ar1  +  0r2  .  (75) 

o 

will  follow  a  normal  distribution  with  variance  o3  where 

=  <x2a\  +  02o*  .  (76) 

If  it  is  recognized  that  a  normal  distribution  is  the  same  as  the 
Maxwellian  distribution  for  a  single  velocity  component.,  then  this  result 
implies  the  distribution  for  the  center-of -mass  velocity  components 
obtained  by  averaging  over  N  molecules  as  m  Eqs.  (68)  -  (70).  The  result 
is  that  this  mean  velocity  follows  a  Maxwellian  velocity  distribution 
appropriate  to  a  "super”  molecule  whose  mass,  ms,  is  given  by 

N  N 

ms  =  <  2W 

i=l  i=l 

(The  temperature  of  the  distribution,  of  course,  is  the  same  as  that  used 
to  select  the  constituent  molecular  velocities.)  Although  Eq.  (77)  is  not 
intuitively  obvious  (to  these  authors,  anyway),  it  does  yield  some  expected 
limits.  If  ail  of  the  weighting  factors  are  the  same,  then  ms  is  the  sum 
of  the  masses  of  the  individual  molecules.  However,  if  one  moieculc's 
weighting  factor  is  much  larger  than  the  others  (resulting  in  the 


^i)2/  (  2Wimi) 


(77] 
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center-of-mass  velocity  of  the  group  being  essentially  equal  to  that 
molecule's  velocity),  then  the  distribution  of  center-of-mass  velocity  is 
the  same  as  the  aistribution  for  that  one  molecule. 

72.3  Molecular  Relative  Velocitv  distribution 

The  relative  velocity  between  an  individual  simulated  molecule 
(referred  to  as  "molecule  j")  of  mass  irtj  with  respect  to  the  center-of-mass 
velocity  of  N  other  molecules  will,  therefore,  have  the  same  distribution 
as  the  relative  velocity  between  that  molecule  and  another  molecule  of  mass 
ras.  It  is  a  well  known  result  that  this  velocity  distribution  is  a 
Maxwellian  distribution  appropriate  to  a  molecule  with  a  reduced  mass,  Mjs, 
given  by 


‘js 


mmg 

m+m,. 


(78) 


Put  in  terms  of  the  chi-square  distribution  which  is  used  extensively  in 

2 

the  Monte  Carlo  model,  u.  (the  square  of  the  relative  velocity  between 

J  S 

molecuie  j  and  the  center-of-mass  velocity  of  the  other  N  molecules),  can 
be  expressed 


(79) 


where  Xtj  is  a  variable  selected  from  a  chi-square  distribution  for  the 
relevant  three  translational  degrees  of  freedom.  (k  is  Boltzmann's 
constant,  ano  T  is  the  (as  yet  undetermined)  temperature.] 

7.2.4  Translational  Energy  of  Relative  Motion 

The  total  translational  energy  of  the  molecules  (which  can  be 
expressed  as  in  Eqs.  (62)  -  (71)  above,  summing  over  all  N+l  molecules)  can 
be  algebraical ly  recast  in  a  form  which  specifically  shows  the  contribution 
of  relative  motion  between  molecule  j  and  the  N  other  molecules. 
Speci ficai ly , 


E 


N+l 


E 


N 


+ 


2 

js 


(80) 
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In  Eq.  (80),  E^  is  the  translational  energy  which  would  result  if  only  the 
N  molecules  were  included  in  the  previous  sums,  and  Ejj+1  is  the  value 
obtained  with  all  of  the  molecules.  The  factor  in  the  difference,  rjj ,  is 
the  "reduced  weighted  mass"  between  molecule  j  and  ail  of  the  other 
molecules,  i.e., 


si+wjmj 


(81  ) 


where  Sj  is  as  defined  in  Eq.  (62),  applying  the  sum  to  the  N  remaining 
molecules.  It  is  crucial  to  note  that  ms ,  as  defined  in  Eq.  (77), 
determines  the  distribution  of  relative  velocities  between  mj  and  the  other 
N  molecules;  but  iq j ,  as  defined  above,  determines  the  amount  of  energy 
associated  with  that  relative  velocity.  Combining  Eqs.  (79)  and  (80)  gives 
the  translational  energy  contribution,  Ef. j ,  as 


E 


tj 


<IkT> 


(82) 


Note  that  this  effectively  gives  a  weighting  factor  associated  with  the 
translational  energy  contribution  of  molecule  j  of  r)j/Mjs- 

7.25  Determination  of  Temperature 

The  internal  energy  associated  with  molecule  j,  Ej,j,  can  be 
represented  simply  by 


EU  =  (ikT,WjXij  '  (83) 

where  Xxj  is  a  variable  selected  from  a  chi-square  distribution  with  the 
number  of  degrees  of  freedom  appropriate  to  molecule  j's  internal  moaes. 
As  discussed  above,  this  process  is  then  repeated  sequentially  for  each 
molecule  in  the  ceil.  Note  that  "N"  in  the  above  relations  refers  to  all 
remaining  molecules  which  have  not  had  their  energies  determined  yet.  This 
means  that  ms,  for  instance,  changes  with  each  molecule  since  it  is  defined 
via  a  sum  over  these  remaining  molecules.  The  last  molecule  has  no 
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transiationai  energy  of  relative  motion  associated  with  it  since  there  are 
no  remaining  molecules  for  it  to  be  moving  with  respect  to.  It  does,  of 
course,  have  internal  energy. 


Summing  ail  of  the  Etj  and  and  equating  them  to  the  known  total 
energy  Etot  (as  given  in  Eq.  (73))  then  determines  the  temperature  of  the 
system.  it  is  noteworthy  that  this  temperature  is  not  determined  by  the 
total  energy  of  the  system,  but  also  by  the  statistical  sampling  process. 
This  is  consistent  with  the  fact  that  any  temperature  could  result  in  the 
particular  observed  velocities;  although  some  temperatures  are  much  more 
likely  to  produce  them  than  others. 

Once  the  temperature  is  defined,  then  the  sequential  relative 

velocities  squared  u^„  (Eq.  (79))  and  Internal  energies  (Eq.  (03))  are 

J  s 

determined.  It  is  then  a  simple  matter  to  go  back  and  apply  these  values 
to  select  individual  molecular  velocity  components  via  the  same  procedure 
described  in  Section  5. 

7.3  The  NumDer  of  Collisions  Required  to  Achieve  Equilibrium 

The  number  of  collisions  required  to  achieve  equilibrium  depends  on 
the  mooei  being  employed  and  the  criterion  for  equilibrium.  (Equilibrium 
is  approached  asymptotically  and,  as  such,  could  be  regarded  as  an  Ideal 
limit  which  is  never  realized.)  The  model  being  employed,  as  discussed  in 
Subsections  2.4  and  5.3  is  that  of  Ref.  4.  In  this  "statistical  collision" 
model,  a  fraction,  a.  of  the  collisions  are  taken  to  be  "perfectly 
inelastic";  that  is,  in  such  collisions  all  translational  and  rotational 
energy  of  the  colliding  molecules  is  made  available  for  distribution  to  the 
post-collision  state  vectors.  taking  into  account  the  number  of 
translational  and  internal  degrees  of  freedom.  The  rest  of  the  collisions 
are  taken  to  be  completely  elastic,  with  no  interchange  taking  place 
between  the  translational  and  rotational  energy  modes.  The  parameter  of 
the  model,  a,  should  be  chosen  to  match  available  data  for  rotational 
relaxation . 

Within  the  context  of  this  model,  the  question  to  be  addressed  is  how 
many  collisions  are  required  in  a  cell  before  the  model  predicts  that  it  is 
essentially  in  equilibrium.  The  question  can  be  made  independent  of  a  it 
it  is  phrased:  "How  many  inelastic  collisions  per  molecule  must  be 
simulated  before  the  cell  can  be  considered  to  be  in  equilibrium?".  This 
quest  ion  is  suitable  tor  direct,  investigation  with  the  model,  and  a  test 
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calculation  was  performed  to  answer  it.  The  test  calculation  indicated 
that  the  equilibration  is  90%  compLet.e  after  approximately  3.08  inelastic 
collisions  (on  the  average)  for  each  molecule.  This  seems  to  represent  a 
reasonable  point  at  which  to  say  that  further  collision  simulation  is 
unnecessary,  although  the  cutoff  is  of  necessity  somewhat  arbitrary.  This 
number  of  inelastic  collisions  serves  as  a  useful  benchmark  m  the 
comparison  of  the  collision  cutoff  and  equilibrium  aftermath  approaches, 
and  it  also  serves  to  define  when  the  application  of  the  equilibrium 
aftermath  approach  is  valid. 

7.4  Method  Comparison 

Test  runs  were  run  wnere  the  collision  cutoff  approach  was  utilized 
for  3.08  inelastic  collisions  per  moiecule.  (Since  a  =  0.2  was  used,  this 
corresponded  to  about  15  total  collisions  per  molecule.)  The  time  required 
to  compute  the  relaxation  via  collisions  was  then  compared  to  the  time 
required  to  utilize  the  equilibrium  aftermath  approach.  The  result  was 
that  the  equilibrium  aftermath  approach  was  almost  an  order  of  magnitude  (a 
factor  of  9)  faster  in  achieving  the  same  result.  This  ratio  will  no  doubt 
vary  with  computer  and  specific  calculation  being  performed;  but  it  is 
nighly  likely  that  the  equilibrium  aftermath  will  always  come  out 
considerably  faster.  It  is  for  this  reason  that  the  method  was  implemented 
in  SOCRATES. 

8.  MOLECULAR  TRANSLATIONS 

As  discussed  in  Section  1,  an  essential  element  of  the  direct 
simulation  Monte  Carlo  method  is  the  periodic  advancement  of  simulated 
molecules  along  their  trajectories.  Formally,  this  is  accomplished  by 
updating  the  position  and  velocity  elements  of  the  state  vector.  The 
specific  procedures  for  doing  this  depend  on  the  coordinate  system  in  which 
the  state  vector  elements  are  represented. 

8.1  Molecular  Translations  in  Cartesian  Coordinates 

In  Cartesian  coordinates,  the  translation  is  very  direct.  Let  x,  y, 
ana  z  represent  the  position  coordinates  and  u,  v,  ana  w  the  corresponding 
velocity  coordinates.  If  initial  and  final  values  of  the  state  vector  are 
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represented  by  a  0  and  1  subscript  respectively,  then  the  updated  state 
vector  elements  corresponding  to  a  translation  through  a  time  step  At  are 


;jj  ver:  by 

x1  =  x0  +  u0At  ,  (84) 

V]  =  Vo  +  v0At  ■  (85) 

■Ix  -  z0  +  w0At  ,  (86) 

Uj  =  U0  (87) 

v,  =  V0  (88) 

ana 

Wj  =  w0  (89) 


8.2  Molecular  Cloning 

When  a  simulated  molecule  is  translated  from  one  cell  to  another,  the 
weignting  factor  for  that  species  will  generally  be  different  in  the  new 
ceil  Since  it  is  the  number  of  real  molecules  rather  than  the  number  of 
simuiateu  molecules  which  must  be  preserved  when  crossing  cell  boundaries 
(statistically,  at  least),  it  is  necessary  to  correct  for  the  distinct 
weighting  factors  (see  Subsection  3.4). 

If  the  weighting  factor  before  translation  is  WQ.  then  the  simulated 
molecule  represents  that  many  real  molecules.  If  the  weighting  factor  in 
the  new  ceil  is  Wj ,  then  Wq/W^  simulated  molecules  would  be  required  to 
represent  the  same  number  of  real  molecules  in  the  new  cell.  If  this  ratio 
were  a  whole  integer,  then  this  could  be  accomplished  by  introducing  that 
many  clones"  of  the  simulated  molecules  in  the  new  cell.  That  is,  Wq/Wj 
simulated  molecules  would  be  placed  in  the  new  cell,  all  with  the  same 
state  vector. 

When  the  number  Wq/Wj^  is  not  an  integer  (the  usual  case,  of  course), 
then  the  cloning  must  be  done  on  a  statistical  basis.  So,  for  instance,  if 
Wq/Wj  were  equal  xo  27,  then  30%  of  the  time  two  clones  would  be  produced, 
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ana  70%  of  rhe  time  three  clones  would  be  proaaced.  Note  that  the  ratio 
may  be  less  than  unity,  ana  the  molecule  may  not  be  introduced  Into  the  new 
ceil  at  ail.  (In  which  case  the  molecule  is  removed  from  the  simulation.) 

Cloning  is  a  necessary  evil  inherent  in  a  system  with  spatially 
varying  weighting  factors.  It  enables  such  a  system  to  maintain  the 
statistically  correct  flux  of  mass  and  momentum  across  cell  boundaries,  but 
it  misrepresents  the  flux  of  randomized  or  thermal  energy.  This  can  be 
seen  by  an  extreme  case  where  a  very  large  number  of  clones  is  produced 
when  a  simulated  molecule  crosses  a  cell  boundary.  The  resulting  molecules 
in  the  new  cell  have  the  correct  mass  and  momentum  flux,  but  since  they  all 
have  precisely  the  same  velocity  they  have  a  null  relative  velocity  arid, 
therefore,  a  zero  temperature.  If  the  weighting  factors  are  not  too 
different  between  adjacent  cells,  then  the  errors  introduced  by  this 
process  are  acceptably  small.  However,  it  does  mean  that  one  cannot, 
arbitrarily  improve  statistics  in  one  portion  of  the  solution  region  by 
selectively  reducing  the  weighting  factors  there.  This  was  a  difficulty 
which  was  encountered  in  the  early  stages  ol  the  direct  simulation  Monte 
Carlo  method  while  trying  to  improve  statistics  along  the  axis  ol 
axtsymmetric  simulations,  since  the  cell  volumes  (and,  therefore,  the 
sample  sizes)  tend  to  be  smallest  on  the  axis. 

As  was  the  case  for  simulated  molecules  produced  via  chemical 
reactions,  it  is  possible  for  the  weighting  factors  between  successive 
ceils  to  be  so  different  that  a  prohibitively  large  number  of  simuiateu 
molecules  would  be  required  to  produce  the  same  number  of  real  molecules. 
Tiie  codes  sense  when  a  disproportionate  number  of  simulated  molecules  are 
being  produced  for  a  given  species  and  cell  and  adjust  the  weighting  factor 
automatically.  As  the  weighting  factor  is  increased,  a  proportionate 
fraction  of  molecules  of  that  species  and  cell  are  removed  from  the 
simulation  in  order  to  keep  the  number  of  real  molecules  properly 
represented.  This  process  enables  the  weighting  factors  to  seek  their  own 
proper  level  without  a  priori  knowledge  of  the  solution.  (Periodically, 
the  ceils  are  examined  to  determine  if  a  certain  species  has  been 
underrepresented  in  terms  of  its  number  of  simulated  molecules.  If  this  is 
found  to  be  the  case,  then  the  weighting  factor  is  decreased,  allowing 
weighting  factors  to  float  downwards  as  well  as  upwards.  It  is  the  danger 
of  weighting  factors  being  too  small,  causing  an  overflow  of  code 
dimensions,  which  is  most  critical,  however.) 
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9. 


INITIAL  AND  BOUNDARY  CONDITIONS 


9.1  Initial  Conditions 

Since  the  direct  simulation  Monte  Carlo  method  is  inherently  an 
unsteady  technique,  an  initial  state  must  be  specified  in  order  to  advance 
the  solution.  For  situations  where  a  steady  state  result  is  desired,  it  is 
obtained  as  the  long-time  solution  to  an  unsteady  problem.  In  this  case 
the  initial  conditions  have  no  effect  on  the  eventual  solution,  but  they 
may  well  have  an  impact  on  the  speed  with  which  that  state  is  achieved. 
For  steady  state  solutions,  SOCRATES  simply  starts  with  an  evacuated 
solution  region.  For  unsteady  solutions,  however,  it  is  necessary  to  start 
with  a  molecular  distribution  which  is  representative  of  the  conditions  at 
the  start  of  the  desired  simulation.  For  SOCRATES  these  conditions 
correspond  to  a  uniform  flow  with  the  translational  and  internal  modes 
being  in  equilibrium.  The  specification  of  the  initial  conditions  for 
unsteady  runs,  therefore,  involves  determining  the  state  vector  elements 
consistent  with  this  condition  for  the  desired  number  of  molecules. 

9.1.1  Number  of  Simulated  Molecules  and  Their  Weighting  Factors 

The  desired  number  of  simulated  molecules  for  each  species  in  each 
cell,  Mc,  .is  an  input  quantity.  (Typically,  simulations  aim  for  a  total 
numoer  ol  molecules  per  cell  in  the  neighborhood  of  20.)  Given  the  initial 
number  uensity  to  be  simulateu  for  a  species,  nj  ,  (which  will  have  been 
automatically  converted  to  internal  dimensions  -  see  Section  3)  the 
weighting  factor  for  species  i  in  a  given  cell  is  simply 


wi 


(90) 


where  V  is  the  ceil  volume.  If  a  species  is  not  initially  present  in  a 
cell,  then  it  is  assigned  an  Initial  weighting  factor  of  zero.  If 
simulatea  molecules  come  into  the  cell,  the  weighting  factor  from  their 
piacc  of  origin  will  bo  used  to  initialize  the  weighting  factor  in  the 
ceil.  As  the  solution  proceeds,  the  weighting  factors  are  automatically 
adjusted  to  keep  the  average  number  of  simulated  molecules  of  each  species 
equal  to  M(.  in  each  cell. 
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9  1 . 2  Initial  Positions 


The  initial  molecules  assigneu  to  a  ceil  should  have  ait  equal 
yrobaoiiity  of  being  piacea  in  any  volume  element  ol  the  ceil.  For  me 
nexaneurai  ceils  of  SOCRATES,  this  simply  involves  selecting  each  of  the 
position  elements  at  random  from  the  range  appropriate  to  tne  ceii  in 
question.  That  is,  the  x  position  is  selected  via  the  equation 
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wnere  x  ana  x  are  the  positions  of  the  x- faces  of  the  cell  in 

m  i  ii  lucix 

question.  The  other  position  elements  are  selected  analogously. 


9.1.3  Initial  Velocity  Components 


The  thermal  velocity  components  for  a  molecule  In  translational 
equilibrium  (neglecting,  for  the  moment,  any  mean  flow  contribution)  should 
be  selected  from  a  normalized  Maxwellian  velocity  distribution.  f'u(v), 
given  by 

f0  ^  a  expi  -(av)2  j/-v¥  ,  (92) 


where 


a  -  -vm/(2R0Tw)  ,  (93) 

m  's  me  species  molecular  weight,  Rq  is  the  universal  gas  constant,  anu  Ttt 
'  s  the  temperature.  Equation  (93)  applies  for  each  of  the  molecular 
velocity  components  and  must  be  sampled  three  times  for  each  molecule  that 
comprises  the  initial  state  of  the  simulation.  A  method  for  directly 
sampling  from  this  distribution  is 

Aj  =  2ir(3  ,  (94) 

AZ  =  ~-V-log(0)  (95) 


v  =  A^sinfAj) 


(96) 
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After  the  thermal  velocity  components  are  determined  for  each  molecule, 
then  any  mean  flow  velocity  is  simply  added  on.  The  velocities  are  then 
transformed  to  internal  units. 

9.1.4  Initial  Internal  Energies 

The  only  remaining  element  of  the  state  vector  to  be  specified  is  the 
internal  energy.  Internal  energies  for  a  gas  in  equilibrium  are 
distributed  according  to  the  normalized  distribution  function  fj  given  by 

f.  =  i-VJj.  ,  (97) 

2«/2r(C/2) 
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represents  the 
in  question,  r 
energy ,  l . e .  , 

2Ej 


number  of  internal  degrees  of  freedom  for  the 
is  the  gamma  function,  and  £  is  a  dimensionless 


(OH) 


where  Ej  is  the  internal  energy.  Equation  (97)  is  a  representation  of  the 
chi  square  distribution  for  C  degrees  of  freedom;  procedures  for  sampling 
from  this  distribution  are  given  in  Appendix  A. 

9.2  Source  Boundary  Conditions 

The  introduction  of  source  molecules  into  the  simulation  is  a  boundary 
condition  which  depends  on  the  specific  model  for  the  source  in  question. 
SOCRATES  includes  a  "core  flow"  source,  which  describes  the  contamination 
that  results  from  the  scattering  of  the  flow  from  a  thruster  back  onto  the 
shuttle.  For  this  source,  it  is  important  to  have  a  description  of  the 
main  exhaust  flow  away  from  a  thruster.  (This  is  to  be  contrasted  with  a 
source  which  describes  the  direct  contamination  of  surfaces  via  impact  of 
exhaust  gases.  Since  the  thrusters  are  not  pointed  directly  at  shuttle 
surfaces,  it  is  the  small  portion  of  the  flow  which  leaves  the  thruster  at 
a  large  angle  from  the  exhaust  centerline  which  is  important  in  this  case.) 

The  plume  gases  expand  upon  leaving  the  exit  plane  and  adopt,  an 
essentially  radial  flow  profile  over  a  distance  which  is  on  the  order  of 


49 


exit  plane  dimensions.  Since  this  distance  is  small  compared  to  the  length 
scales  of  the  interaction  of  the  plume  with  the  atmosphere,  it  is 
appropriate  to  replace  the  nozzle  by  a  point  source  of  exhaust  molecules 
traveling  at  their  thermoaynamic  limiting  speed  with  an  undisturbed  number 
density  distribution  given  by 


(99) 


where  B  gives  the  axial  number  density  decay  and  e  represents  the  angie 
from  the  thrust  axis.  The  r  appearing  in  Eq.  (99)  is  the  spherical  radius, 
giving  the  total  distance  from  the  source.  The  particular  form  rur  f(e) 
that  is  used  in  SOCRATES  is  an  asymptotic  form  of  that  proposed  by  Btook.7 
namely 

1(0)  =  exp{-X2(l  -  cos(e)]}  ,  (100) 
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and 


1  *  ,2U« 
— n„ApX6 — 
2ir  e  e  um 
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in  the  above  relations,  ue ,  Mc,,  np  and  Ae  denote  the  exit  velocity, 
Mach  number,  number  density  ana  area,  respectively;  en  is  the  nozzle 

7.  Btook,  J.  W. ,  "Far  Field  Approximat ion  for  a  Nozzle  Exhausting  into  a 
Vacuum",  Journal  of  Spacecraft  and  Rockets,  626  (1969). 
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divergence  half  angle,  and  um  and  T  are  the  thermodynamic  limiting  speed 
ana  the  ratio  of  specific  heats. 


in  a  solution  time  step  £tm.  NelleAeAtm  real  molecules  are  introduced. 
Eacn  moiecuie  is  assigned  an  angle  e  which  is  chosen  to  be  consistent  with 
r.q.  (100)  via 


e  =  cos  1 t 1  +  C1log(l  -  C2#)] 

(105) 

wnere 

1 

(106) 

and 

C2  =  1  -  exp(-2\2) 

(107) 

An  azimuthal  angle  is  selected  at  random  for  the  moiecuie,  and  then 
the  resulting  velocity  is  represented  in  the  basic  Cartesian  coordinates 
utilized  by  SOCRATES.  The  moiecuie  is  then  advanced  from  the  source 
appropriate  to  a  speed  ol  um  and  a  time  increment  which  is  a  random 
fraction  of  dtm.  The  process  is  repeated  until  the  proper  number  or 
simulated  moiecuies  of  each  exhaust  species  have  been  introduced. 

9.3  Atmospneric  3oundary  Condition 

The  atmospheric  boundary  condition  for  SOCRATES  is  that  molecules 
snouid  be  introduced  into  the  solution  region  from  the  outer  boundaries  in 
sucn  a  way  as  to  simulate  the  undisturbed  ambient  flow  outside  of  the 
solution  region. 

9.3.1  Molecular  Flux  Across  a  Surface  Element 

The  relations  for  molecular  flux  across  an  infinitesimal  surface 
element  are  given  in  Ref.  1.  If  q  is  the  molecular  flux  (molecules  per 
unit  area  per  unit  time)  crossing  a  given  surface  element,  then  q  is  given 
by 

n  _ 

<1  =  — <exp(- +  w(l+erf(w)j}  ,  (108) 

Ci\ 
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woere 


A  =  -Vn/2  R0Tb 


(109) 


arid 


w  =  Aumcos(e)  (HO) 

In  these  relations,  n^  ann  T  represent  the  ambient  number  density  anu 
temperature  for  the  species  in  question;  m  represents  its  molecular  weight; 
uo>  represents  the  mean  flow  velocity;  and  8  is  the  angle  between  the  inward 
surface  normal  and  the  mean  flow  direction.  The  flux  given  by  Eq.  (100)  is 
non-zero  for  all  values  of  e,  reflecting  the  distribution  of  molecular 
velocities.  However,  it  does  become  exponentially  small  for  large  negative 
w.  (Note  that  these  relations  must  be  applied  on  a  spec ies-by- species 
basis;  each  species  has  a  different  spread  in  its  velocity  distribution  by 
virtue  of  its  different  molecular  weight.) 

The  application  of  Eq.  (108)  to  the  flat  surfaces  comprising  the 
SOCRATES  outer  boundary  is  direct,  since  e  does  not  change  along  the  face. 
The  total  number  of  molecules  to  introduce  for  a  time  step  of  Atm  is  simply 
qdtmA,  where  A  is  the  area  of  the  flat  face  in  question.  Since  the  flux  is 
constant  over  a  flat  face,  each  position  on  the  face  is  equally  likely  as  a 
point  for  molecular  entry.  Hence,  for  a  flat  face,  the  starting  molecular- 
position  can  be  simply  obtained  by  selecting  a  point,  at  random  on  the  face. 

9.3.2  Incoming  Molecular  Velocity  Components 

For  each  molecule  that  is  introduced,  a  local  orthogonal  coordinate 
system  is  set  up  such  that  one  direction  is  in  the  direction  of  the  inward 
surface  normal.  Velocity  components  are  first  determined  In  terms  of  this 
local  coordinate  system  and  then  transformed  to  the  main  code  coordinate 
system.  In  the  local  coordinate  system,  the  velocity  components  parallel 
to  the  surface  are  determined  as  discussed  above  for  molecules  in  the 
initial  condition.  The  inward  component  of  velocity  must  be  selected  in 
proportion  to  the  distribution  h(v)  given  by 

h(v)  =  v  exp{-la(v  -  <v>)]2)  ,  (in) 
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wuere  <v>  is  the  component  of  the  mean  flow  velocity  in  the  inward  normal 
oirecrion  ana  a  is  as  given  in  Eq.  (03).  It  is  possible  for  <v>  to  be 
negative,  but  all  incoming  molecules  must  have  a  positive  v  value  by 
definition.  Hence,  this  distribution  is  only  sampled  for  positive  values 
of  v.  The  sampling  is  done  via  the  acceptance-rejection  technique. 


10.  STATISTICAL  SAMPLING  OF  OUTPUT 
10.1  General  Considerations 

It  is  safe  to  say  that  the  molecular  state  vectors  as  they  exist  in 
the  computer  do  not  comprise  the  usual  desired  output  of  the  procedure. 
With  rare  exceptions,  it  is  the  macroscopic  quantities  such  as  temperature, 
density,  mean  flow  velocity,  etc.  which  are  of  interest  -  not  the 
microscopic  quantities  represented  by  the  state  vector  of  an  individual 
simulated  molecule.  The  generation  of  the  desired  output  requires  that  the 
macroscopic  quantities  of  interest  be  represented  in  terms  of  statistical 
sums  of  r.he  available  microscopic  quantities;  and  it  is  the  main  purpose  of 
this  section  to  present  these  correspondences.  All  sums  are  kept  in  terms 
oi  "real"  molecules  and  events,  i.e.,  the  current  weighting  factors  are 
induced  in  the  sums.  This  is  essential  since  the  weighting  factor 
determines  the  statistical  importance  of  a  given  molecule.  Since  the 
weighting  factors  are  dynamically  and  unpredictably  adjusted  as  the 
solution  progresses,  it  would  not  be  possible  to  go  back  and  add  in  the 
effect  of  weighting  factors  a  posteriori. 

In  general,  it  must  be  decided  ahead  of  time  exactly  what  output  is 
ucsireu  from  the  code,  and,  therefore,  what  statistical  sums  should  be  kept 
to  generate  it.  There  is  a  vast  amount  of  potential  information  in  the 
simulation,  and  it  is  not  reasonable  to  store  all  possibly  interesting 
quantities  in  all  runs.  On  the  other  hand,  it  is  wasteful  to  completely 
rerun  a  case  just  because  the  user  decides  there  was  an  additional  quantity 
he  was  interested  in.  The  selection  of  output  for  a  given  run,  therefore, 
unavoidably  requires  user  judgment.  Once  the  user  has  decided  upon  the 
required  output,  the  determination  of  which  statistical  sums  are  required 
is  done  automatically  by  the  code.  Care  is  taken  to  make  sure  that  a 
statistical  sum  is  not  duplicated  internal ly  if  it  is  required  by  more  than 
one  requested  output  quantity. 
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Some  initial  words  of  caution  are  required.  By  its  nature,  the  direct 
simulation  Monte  Carlo  method  works  with  far  fewer  molecules  than  nature 
(toes,  and  it,  therefore,  exhibits  consioerably  greater  statistical 
variation  in  its  macroscopic  predictions.  To  reduce  these  variations,  the 
coae  is  run  repeatedly  for  the  same  case,  increasing  the  statistical  base 
from  which  the  macroscopic  output  is  derived.  Useful  results  can  usually 
be  obtained  with  a  modest  computational  effort.  However,  this  statement 
must  be  tempered  by  a  realization  of  the  convergence  rate  for  Monte  Cano 
sampling.  Basically,  the  statistical  error  in  the  output  converges  as  one 
over  the  square  root  of  the  sample  size  (or  run  time).  Hence,  if  a 
solution  looks  good,  but.  the  user  decides  he  would  like  one  more 
significant  digit  (i.e.,  he  would  like  trie  statistical  error  to  be  reduced 
to  0.1  times  its  current  value)  it  would  require  that  the  run  time  be 
increased  by  a  factor  of  100!  It  can  be  seen  that  the  desire  for  more 
accuracy  can  quickly  turn  the  most  efficient  code  into  a  money  gobbling 
nightmare.  When  using  a  Monte  Carlo  technique,  one  must  acceyi  some 
statistical  scatter  in  the  output. 

10.2  Sampling  of  Instantaneous  Volumetric  Output  Quantities 

Instantaneous  volumetric  output  quantities,  such  us  density, 
temperature  and  velocity,  can  be  determined  by  examining  the  molecular 
stale  vectors  at  a  particular  time  in  the  simulation.  The  code  pauses  m 
the  simulation  and  uses  the  molecular  state  vector  elements  to  aud  values 
to  the  statistical  sums  appropriate  to  the  various  cells  and  the  particular 
time  ihat  it.  paused.  It.  then  proceeds  with  the  simulation  until  the  next 
sampling  time.  As  the  code  goes  through  its  successive  runs,  it  stops  at 
the  same  points  in  the  simulation  every  time  and  adds  to  the  statistical 
base  tor  the  sums.  (For  steady  state  cases,  it  simply  does  it  repeatedly 
after  the  initial  transient  has  died  down.)  The  items  listed  below,  witii 
their  statistical  definitions,  are  selectable  as  output  requests  in 
SOCRATES.  Summations  are  performed  over  all  applicable  simulated 
molecules,  which  include  Nruri  separate  runs. 

•  TOTAL  NUMBER  DENSITY 


n 


VN 


run 


•  MEAN  MOLECULAR  WEIGHT 


(112) 
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•  X  VELOCITY  COMPONENT 


SS 


•  y  VELOCITY  COMPONENT 
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•  i  VELOCITY  COMPONENT 


•  OVERALL  TRANSLATIONAL  TEMPERATURE 
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•  TRANSLATIONAL  TEMPERATURE  IN  jTH  DIRECTION 


1  S8 
-( s7  -  — ) 


RqSj 


>6 


•  INTERNAL  MODE  TEMPERATURE 
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wnere  the  indicated  sums,  S„,  are  defined  by: 
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S10  =  y  •  (129) 

4mm*. 

i 

With  the  exception  of  Eq.  (113),  all  of  the  above  quantities  can  also  be 
defined  and  calculated  for  any  specified  species.  The  sums  are  the  same 
except  that  only  molecules  of  that  species  are  considered.  Before  printing 
output  quantities,  they  are  always  transformed  to  standard  dimensions  from 
the  internal  dimensionless  variables. 
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10.3  Sampling  of  Time  Averaged  Output  Quantities 


Some  additional  quantities  of  interest  are  not  sampled  at  a  separate 
sampling  time  as  described  above,  but  rather  as  the  simulation  evolves. 
Examples  of  such  quantities  are  collision  rates,  reaction  rates,  mean 
velocities  between  molecules,  etc.  For  the  most  part,  tnese  quantities 
depend  on  the  relative  state  of  more  than  one  type  of  molecule,  ana  they 
are  by  their  nature  expressed  as  average  values  over  a  finite  time 
interval .  The  formulas  for  calculating  these  quantities  are  no  more  than 
event  counters,  ana  will  not  be  included  here.  The  following  quantities 
are  currently  available  as  output: 

•  Mean  Relative  Velocity  Between  any  Two  Species; 

•  R.M.S.  Deviation  of  Mean  Relative  Velocity  Between  any  Two 
Species ; 

•  Mean  Product  ol  Cross  Section  Times  Relative  Velocity  Between 
any  Two  Species; 

•  Collision  Rate  Between  any  Two  Species; 

•  Reaction  Rate  for  any  Chemical  Reaction; 

•  Reaction  Rate  for  any  Photochemical  Reaction; 

•  Flux  Rate  for  any  Species  on  any  Surface  Element. 

The  sampling  for  ail  but  the  last  of  these  quantities  occurs  in  the 
collision  simulation  routines.  As  pairs  are  considered  as  possible 
collision  partners,  statistics  are  kept  to  generate  the  first  three 
quantities.  Statistics  on  collisions  and  reactions  are  kept  as  they  occur, 
ana  the  last  quantity  is  determined  in  the  molecule  advancement  routines. 


11.  SURFACE  DEFINITIONS  AND  INTERACTIONS 

An  essential  element  of  the  contamination  problem  is  the  presence  ol 
so  no  surfaces  in  the  flow  field.  This  section  discusses  how  the  surfaces 
are  represented  ana  how  the  interactions  of  the  gas  molecules  with  the 
surfaces  are  recognized,  simulated,  and  recorded. 
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11.1  Snuttle  Representation 


A  simplified  version  of  the  shuttle  has  been  constructed.  This  model 
was  intentionally  abbreviated  for  the  initial  computations.  The  model 
description  is  simply  a  data  file,  which  is  easily  replaced  by  available, 
more  detailed,  shuttle  models.  The  initial  model  for  the  shuttle  geometry 
was  designed  to  form  a  completely  closed  (i.e.,  no  "holes"), 
non-overlapping  surface  which  approximates  the  shuttle  geometry  with  a 
minimum  number  of  surface  elements.  The  surface  elements  are  simple 
geometric  shapes  such  as  rectangles,  triangles,  disks,  cylinders,  and 
cones.  This  first  model  employs  four  surface  types  with  a  total  of  11 
surface  elements.  In  particular,  the  wings  are  represented  by  triangular 
planes  which  currently  have  no  thickness,  but  necessarily  have  a  top  and 
bottom.  The  tail  is  modeled  using  a  combination  of  four  triangular  planes, 
the  shuttle  body  as  the  outer  surface  of  a  cylinder,  the  shuttle  nose  as  a 
cone,  ana  the  aft  end  of  the  shuttle  as  a  disk.  The  model  is  specified  in 
cartesian  coordinates  with  the  origin  placed  along  the  axis  of  the  cylinder 
ac  the  center  of  the  shuttle.  This  preliminary  model  is  shown  in  front, 
top,  anu  side  views  in  Figures  4-6,  respectively. 
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SIDE  VIEW 
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gure  6.  A  Side  View  of  the  Crude  Shuttie  Mode]  Designed  for 
Testing  of  the  SOCRATES  Modei. 
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11.2  Determination  of  Surface  Intersections 


The  interaction  of  species  with  the  shuttle  is  a  crucial  portion  ol 
the  contamination  model,  and  it.  has  two  distinct  facets: 

1)  Calculating  the  point  in  space  ana  time  at  which  a 
contaminant  molecule  makes  contact  with  a  shuttle  surface: 

2)  Characterizing  what  happens  to  the  molecule  after  contact 
(e.g. ,  adsorption,  specular  reflection,  diffuse  reflection, 
etc. ) . 

This  subsection  deals  with  the  development  of  an  algorithm  for  the 
first  point  above.  The  calculation  of  an  intersection  point,  while 
conceptually  straightforward,  is  a  potential  source  of  considerable 
computational  effort.  SOCRATES  can  calculate  the  intersection  point  in 
space  and  time  for  a  molecule  starting  from  an  arbitrary  position  and 

velocity  for  each  of  the  simple  geometric  shapes  used  in  the  shuttle 
description.  The  intersection  routines  also  return  the  local  triple  of 
unit  vectors  at  the  intersection  point  which  is  useful  for  the  calculation 
of  surface  reflections.  The  procedure  will  be  illustrated  for  the  case  of 
a  rectangular  surface  element.  The  surface,  as  shown  schematically  in 
Figure  7,  is  defined  by  the  following  quantities: 

a)  A  vector,  rs,  giving  the  absolute  location  of  the  "key 

vertex"  of  the  rectangle  in  code  coordinates: 

b)  An  orthonormal  triple  of  unit  vectors  which  define  the 

orientation  of  the  surface.  ij  and  i2  define  the  directions 
from  the  key  vertex  to  the  two  adjacent  vertices  of  the 
rectangle  and  i3  is  the  outward  surface  normal.  A  right 
handed  coordinate  system  is  used,  so 

i3  =  i2  x  i2  ;  (130) 

c)  The  lengths,  and  12,  of  the  two  sides.  (See  Figure  7.) 

If  a  molecule  has  a  position,  rm,  and  a  velocity,  vm,  then  the 

analysis  for  intersection  proceeds  as  follows: 

1)  The  component  in  the  i3  direction  of  the  moiecuie's  position 
ana  velocity  relative  to  the  key  vertex,  x3  and  v3,  are 
computed  via 

x3  =  V<rm-?s>  (131) 
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and 


v3  =  i3.Cm  (132) 

2)  In  order  tor  an  intersection  to  take  place  on  the  proper  side 
of  the  rectangle,  x3  must  be  positive  and  v3  must  be 
negative.  If  these  criteria  are  not  met,  no  further  analysis 
is  performed. 

3)  If  the  above  criteria  are  met.  the  intersection  with  the 
plane  of  the  rectangle  takes  place  at  a  time  increment,  At^, 
given  by 

x3 

Atf  - -  (133) 

1  v3 

4)  The  position  ol  the  intersection  point,  Xjj  and  ,  relative 
to  the  key  vertex  is  then  given  by 


Figure  7.  An  illustration  of  the  Quantities  Used  to  Calculate  a 
Molecular  Intersection  with  a  Rectangular  Plate. 
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(135) 


x2i  =  +  -  rs) 

5)  An  intersection  with  the  rectangle  occurs  if  anti  only  if  (0  £ 
xli  ^  i1)  and  (0  £  x21  £  «2). 

The  procedure  for  the  othe^  surface  types  is  not  given  here,  hut  it  is 
quite  similar.  Each  surface  is  defined  by  a  location  vector,  a  triple  of 
orientation  vectors,  and  a  few  pieces  of  auxiliary  information  which  are 
specific  to  the  surface  type.  The  use  of  simple  geometric  shapes  allows 
the  checks  for  intersection  for  all  of  the  surfaces  to  be  maae 

expeditiously. 

There  wiil  he  thousands  of  molecules  in  a  simulation,  and  each  one  of 
these  molecules  is  advanced  along  its  trajectory  at  every  solution  time 
stop.  Some  check  for  intersection  must  be  made  for  each  molecule  at  every 
time  step.  The  shuttle  is  modeled  as  a  combination  of  several  simple 

surfaces.  Although  the  initial  model  does  not  involve  a  large  number  of 
surfaces,  it  is  an  obvious  growth  path  for  the  contamination  mooel  to  use  a 
more  ana  more  sophisticated  model  of  the  shuttle  itself.  (Other  models9 

have  used  hundreds  of  distinct  surfaces  to  describe  the  shuttle.)  Each 

molecule  may  be  checked  for  possible  intersection  with  every  surface 
element  at  every  step.  It  is  not  even  valid  to  stop  checking  for 
intersections  when  one  is  found,  since  it  is  the  first  intersection  point 
along  the  molecular  path  that,  is  the  one  of  interest;  there  may  be  more 
than  one  mathematical  intersection.  Hence,  it  is  desirable  to  have  an 
algorithm  which  does  not  suffer  greatly  from  a  large  number  of  surfaces. 

A  concept  to  speed  up  the  calculation  of  surface  intersections  was 
implemented.  An  element  was  added  to  the  state  vector  to  indicate  the  time 
at  which  a  given  molecule  will  experience  a  surface  collision  if  its 
current  trajectory  is  not  altered.  The  element  is  used  as  follows: 

1)  Whenever  a  molecule  is  introduced  into  the  simulation,  this 
eiement  is  set  to  zero.  This  serves  as  a  flag  indicating 
that  a  possible  surface  intersection  has  not  yet  been 
calculated  for  this  molecule. 


8.  Hetrick,  M.  and  strange- Jensen ,  D.,  "Shuttle  Computer  Model  for  On- 
Orbit  Contamination  Analysis",  Proceedings  of  the  10th  JANNAF  PI  time 
Technology  Meeting,  CPIA  Publication  291,  175  (1977). 
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2 )  During  collision  sampling,  whenever  a  molecule  has  its 
trajectory  changed,  the  state  vector  element  for  surface 
intersection  is  reset  to  zero.  This  flags  the  molecule  to 
have  its  possible  surface  intersection  recomputed  when  it  is 
advanced  along  its  trajectory. 

3)  The  routines  which  advance  a  molecule  along  its  trajectory 
examine  this  element.  If  it  is  zero,  then  all  surfaces  are 
checked  for  possible  intersection.  If  an  intersection  is 
found,  then  trie  time  at  which  the  intersection  will  take 
place  is  put  in  the  state  vector  element.  If  it  is 
determined  that  the  current  trajectory  will  not  intersect  any 
surfaces,  then  the  value  of  10^°  (a  computer  approximation 
for  infinity)  is  put  in  the  element. 

4)  If  it  is  known  that  a  molecule  will  not  intersect  a  surface 
within  the  time  interval  corresponding  to  the  moleculiir 
advancement,  then  the  molecule  is  simply  moved  along  its 
trajectory  without  further  checking  of  surfaces.  This  will 
be  the  case  for  the  large  majority  of  molecules  which  are 
inspected . 

5)  If  a  molecule  does  intersect  a  surface  within  the  current 

time  interval,  then  it  is  advanced  to  the  point  of 
intersection.  The  stare  vector  corresponding  to  the 

post-ref lection  conditions  are  calculated,  and  the  element 
corresponding  to  surface  collision  time  is  reset  to  zero. 

The  molecule  is  then  advanced  along  its  new  trajectory  for 
the  remainder  of  the  time  interval,  allowing  for  any  new 
reflections  which  may  occur. 

Another  concept  that  has  been  developed  for  speeding  up  the 
calculation  of  surface  intersections  is  to  surround  many  surface  elements 
with  an  artificial  surface  such  as  a  sphere.  If  a  molecule  starts  on  the 
outside  of  the  sphere  and  doesn't  penetrate  it,  then  it  cannot  hit  any  of 
the  surface  elements  within  the  sphere.  In  this  manner,  the  calculation  of 
intersection  for  many  surface  elements  can  frequently  be  replaced  by  the 
cajcuiatlon  of  one  intersection.  (If  the  sphere  is  penetrated,  of  course, 
then  the  detailed  calculations  must  then  be  carried  out.  The  expectation 
is,  however,  that  a  large  fraction  of  molecules  will  not  need  the  detailed 
analysis.)  This  concept  will  be  implemented  in  future  code  versions  if 
computationally  required. 
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11.3  Surface  Reflections 


Routines  were  written  to  describe  a  diffuse  reflection  of  a  moiecuie 
from  a  surface  after  complete  accommodation.  This  is  felt  to  be  the  most 
reasonable  physical  model,  so  it  is  the  natural  initial  choice.  Other 
options  will  be  added  as  the  model  is  expanded. 

11.4  Surface  Statistics 

The  purpose  of  a  contamination  code  cannot  be  served  unless  the  llux 
to,  and  buildup  of,  contaminants  on  the  surface  elements  can  be  described. 
The  ability  to  keep  statistics  for  species  fluxes  on  the  various  surface 
elements  Is  part  of  the  code.  The  requests  are  input  with  the  surface 
element  definitions  and  involve  listing  the  species  for  which  surface  flux 
information  is  desired  for  the  given  surface  element.  The  code 
automatical ly  sets  up  the  required  storage  location,  advances  a  counter 
whenever  a  molecule  of  a  selected  species  strikes  the  surface  element,  and 
generates  output  to  show  the  derived  species  flux.  As  part  of  this 
calculation,  it  was  necessary  to  calculate  the  areas  of  each  surface 
element,  so  the  derived  flux  can  be  given  as  a  number  of  molecules  per  unit 
time  per  unit  area. 

11.5  Interface  of  Shuttle  Model  with  Calculational  Grid 

The  melding  of  the  three-dimensional  Cartesian  cell  grid  structure 
with  the  arbitrary  geometry  of  the  space  shuttle  orbiter  (or,  more 
generally,  any  spacecraft)  posed  a  problem.  The  difficulty  arose  from  the 
fact  that  the  cell  structure  Is  not  molded  around  the  orbiter,  so  the 
boundaries  of  the  orbiter  do  not  correspond  to  cell  boundaries. 

The  cells  are  used  for  two  purposes  in  the  simulation:  1)  to  provide 
positions  for  flow  field  output  quantities  and  2)  to  define  the  location  of 
possible  collision  partners  for  moiecuies.  For  either  purpose,  but 
especially  for  the  second  one,  it  is  desirable  that  macroscopic  properties 
ao  not  change  appreciably  across  the  cell.  This  is  because  the  only 
spatial  requirement  on  two  molecules  to  qualify  as  collision  partners  is 
that  they  lie  within  the  same  cell;  if  the  cell  is  uniform,  then  it  is 
argued  that  a  sample  moiecuie  could  equally  well  be  found  anywhere  within 
the  ceil,  and  its  precise  location  in  the  cell  is  ignored  in  the  collision 
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sampling  process.  (This  assumption  is  not  made  arbitrarily,  of  course;  it 
results  in  a  substantial  computational  simplification.  See  Section  6  for 
oetai is .  ) 

There  is  a  problem  with  this  assumption  when  pieces  of  the  shuttle 
penetrate  into  a  cell.  The  general  collision  sampling  procedure  would 
allow,  for  instance,  a  molecule  above  a  wing  to  collide  with  a  molecule 
oeiow  a  wing.  It  is  an  inherent  contradiction  to  assume  that  the  contents 
of  a  given  ceil  are  uniform  and  that  a  piece  of  the  shuttle  penetrates  the 
cell  since  the  shuttle  piece  defines  a  length  scaie  on  the  order  of  ceil 
uimensions . 

The  resolution  that  was  achieved  was  to  analyze  the  cell  arid  shuttle 
geometries  to  tag  those  cells  which  contained  pieces  of  the  shuttle  and  to 
disallow  collisions  in  those  cells.  As  the  ceils  become  smaller,  the 
neglected  collisions  become  insignificant,  so  this  is  formally  a  source  of 
error  associated  with  finite  grid  size  (an  inherent  part  of  any  such 
calculation).  It  was  judged  better  to  not  allow  the  very  small  number  of 
legitimate  collisions  in  these  cells  than  to  allow  momentum  transfer 
between  molecules  separated  by  a  solid  surface.  It  is  important  to 
emphasize  that  this  approximation  only  has  to  do  with  molecular  collisions. 
Direct  contamination  from  a  source  to  a  surface  still  occurs,  since  it  has 
nothing  to  qo  with  the  cell  structure  at  all. 

11.6  Baclc-To-BacK  Surfaces 

A  shuttle  surface  such  as  a  wing  is  considered  to  have  negligible 
thickness  in  the  first  shuttle  model,  so  it  is  modeled  as  back-to-back 
pianar  segments;  one  for  each  of  the  two  outward  normal  directions.  Due  to 
the  neglect  of  surface  thickness,  when  a  molecule  intersects  with  such  a 
back-to-back  surface,  it  can  potentially  bounce  back  and  forth  between  the 
two  surfaces  indefinitely  without  ever  moving.  The  code  recognizes  this 
situation  and  disallows  it. 

12.  SAMPLE  CALCULATIONS 
12.1  Case  Descriptions 

In  order  to  check  out  the  code  and  demonstrate  some  of  its  current 
capabilities,  sample  calculations  were  undertaken.  The  calculations  were 
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tor  cases  with  the  simplified  shuttle  model  flying  at  altitudes  of  200. 
250,  ana  300  kilometers  at  a  velocity  of  7  km/s.  The  shuttle  was  taken  to 
be  flying  at  a  normal  aircraft  orientation  (i.e.,  nose  into  the  oncoming 
stream)  and  firing  an  HCS  thruster  upward  from  its  nose.  The  geometry  of 
the  calculations,  with  the  coordinate  system,  is  depicted  in  Figure  8. 
These  coordinates  are  important  for  the  understanding  of  subsequent 
resui ts . 


M 


Figure  tt.  A  Schematic  of  the  Sample  Shuttle  Problem  Showing  the 
Coordinate  System  and  Orientation  of  the  Calculation. 
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The  calculation  were  performed  for  a  steadily  firing  engine  of  860 
pounds  thrust,  with  an  exhaust  composition  as  shown  in  Table  2.  This 

composition  was  chosen  to  show  the  specific  effect  of  exhaust  species 

molecular  weight  on  the  development  of  the  contamination  cloud.  Hence, 
light  (H2),  medium  (H20)  and  heavy  (C02)  molecules  were  explicitly  carried 
in  the  simulation,  and  all  other  species  were  grouped  together  into  a 

species  called  "Other",  with  a  molecular  weight  given  by  the  mean  of  all 

the  remaining  exhaust  species. 


Table  2.  Exhaust  Composition  for  Sample  Calculations. 


Species 

Mol .  Weight 

Mole  Fraction 

h2 

2.00 

0.1800 

h2o 

18.00 

0.3284 

co2 

44.00 

0.0472 

Other 

27.35 

0.4444 

The  atmospheric  and  calculational  parameters  used  for  the  three 
altitudes  are  shown  in  Table  3.  The  mean  free  path  is  actually  a  code 
output  variable,  but  it  is  an  important  reference  quantity  so  it  is 
included  in  the  table.  It  should  be  noted,  however,  that  this  mean  free 
path  is  for  the  undisturbed  atmosphere,  and  the  plume-atmosphere 
interaction  region  is  characterized  by  a  smaller  value. 

An  artifice  was  used  to  get  better  statistics  on  the  atmospherically 
scattered  molecules.  An  unscattered  and  scattered  version  of  each  species 
(H2(U),  H2(S),  etc.)  was  defined,  and  gas  kinetic  reactions  of  the  form: 

H2(IJ)  +  M  -*  H2(S)  +  M 

H20((J)  +  M  -»  H20(S)  +  M 

etc.  were  used.  Since  the  contaminant  source  introduced  unscattered 
molecules,  the  flow  field  level  of  H2(U),  for  instance,  then  corresponded 
to  the  portion  of  the  molecular  hydrogen  which  had  not  yet  experienced  a 
collision,  since  any  collision  would  result  in  the  H2(U)  being  transformed 
into  H2(S).  Similarly,  the  H2(S)  represented  the  density  of  scattered 
hydrogen,  and  the  total  hydrogen  density  was  simply  the  sum  of  H2(U)  and 
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Tablo  3.  Atmospheric  and  Calcuiational  Parameters  for 
Sample  Calculations. 


Quantity 

200 

Altitude 

250 

(km) 

300 

Total  Density  (molecuies/cm3 ) 

7.2X109 

1.9x10s 

6.5x10s 

0  Mole  Fraction 

0.58 

0.74 

0.85 

N2  Mole  Fraction 

0.42 

0.26 

0.15 

Temperature  (K) 

855. 

941  . 

976. 

Mean  Free 

Path  (m) 

325. 

1410. 

4540. 

Minimum  X 

(m) 

-500. 

-3000. 

-5000. 

Maximum  X 

(m) 

2000. 

6000. 

10000. 

Minimum  Y 

(m) 

-1000. 

-3000. 

-5000. 

Maximum  Y 

(m) 

1000. 

3000. 

5000. 

Minimum  Z 

( m ) 

-500. 

-3000. 

-5000. 

Maximum  Z 

(m) 

2000. 

6000. 

10000. 

Number  of 

Cells 

2700. 

3888. 

3888. 

H2(S).  The  advantage  ot  this  artifice  is  that  it  allows  the  contamination 
from  scattered  species  to  be  calculated,  even  though  the  unscattereu 
species  densities  are  much  greater  than  the  scattered  species  densities  in 
regions  adjacent  to  the  firing  engine.  By  separating  out  the  important 
scattered  portion,  which  would  otherwise  be  dominated  by  the  unscattered 
portion,  better  statistics  are  obtainable  on  the  desired  quantities. 

12.2  Contamination  Cloud  Results  at  200  Kilometers 

The  presentation  of  quantities  as  a  function  of  three  spatial 
dimensions  is  always  somewhat  difficult.  The  approach  that  has  been  taken 
is  to  show  isodensity  contours  for  various  planes.  The  coordinates  defined 
in  Figure  8  are  used,  and  emphasis  is  placed  on  the  scattered  species, 
since  they  are  the  contamination  source  which  is  being  accurately 
represented  in  this  calculation. 
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Figure  9  shows  isoaensity  contours  tor  the  scattered  H2  molecules  for 
a  piane  perpendicular  to  the  oncoming  wind  direction  at  the  shuttle 
location.  The  molecules  are  more  numerous  above  the  shuttle,  as  is  to  be 
expected  from  the  direction  of  the  firing,  but  a  significant  density  has 
made  it  beneath  the  shuttle.  Figure  10  shows  the  same  quantity  for  a  plane 
0.5  kilometers  in  front  (upstream)  of  the  shuttle.  The  contours  have 
spread  out,  indicating  a  more  even  distribution.  Figure  11  shows  the 
situation  at  1.5  km  upstream,  where  the  scattered  H2  density  has  become 
essentially  constant.  The  contours  are  ragged  since  small  statistical 
fluctuations  become  very  important  in  defining  contour  locations  for  a 
nearly  constant  function. 

Figures  12  and  13  show  isodensity  contours  of  scattered  H2  and  C02, 
respectively,  for  a  plane  at  a  constant  Z  location  of  1  km  above  the 
shuttle.  In  these  figures  the  wind  is  approaching  from  the  right.  The 
comparison  of  the  two  figures  reveals  more  widely  spaced  contours  for  the 
scattered  H2  than  for  the  heavier  C02.  These  figures  quantitatively 
demonstrate  what  could  have  been  predicted  qualitatively;  namely,  that  H2 
is  more  effective  at  traveling  upstream  than  the  other  two  molecules.  This 
is  partially  because  H2  has  a  smaller  cross  section  (and,  therefore,  a 
greater  mean  free  path)  than  the  other  molecules,  but  mainly  because  its 
smaiier  molecular  weight  results  in  greater  molecular  velocities.  It  makes 
xhe  most  of  the  time  it  has  between  collisions.  This  point  is  made  more 
directly  in  Figure  14,  where  the  number  densities  of  the  three  scattered 
species  are  shown  at.  a  location  of  1  km  above  the  shuttle  as  a  function  of 
the  upstream  coordinate.  Within  the  statistical  scatter,  there  is  no 
apparent  difference  in  the  upstream  decay  rate  for  H20  and  C02,  but 
relative  to  these  molecules,  H2  gains  approximately  an  order  of  magnitude 
in  density  at  an  upstream  distance  of  2  kilometers.  (At  this  distance  it 
is  probably  comparable  to  the  atmospheric  concentration  of  H2 ,  although  no 
atmospheric  H2  was  included  In  the  calculation.  The  point  that  is  being 
made,  however,  is  simply  that  light  molecules  are  better  at  making  it 
upstream  and  then  being  blown  back  into  the  shuttle  area.) 

Some  features  of  the  velocity  distribution  function  for  scattered  H2 
molecules  is  shown  in  Figures  15  and  16.  Figure  15  displays  the  mean 
upward  (Z)  velocity  component  for  the  plane  at  X=0.  This  velocity  is  in 
the  neighborhood  of  5  km/s  at  a  distance  of  2  km  above  the  shuttle,  but. 
becomes  negative  at  the  shuttle.  The  spread  in  the  scattered  H2  velocity 
distribution  is  characterized  by  the  translational  temperature  shown  in 
figure  16.  Near  the  shuttle  the  temperature  is  relatively  low  (considering 
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Figure  9.  A  Contour  Plot  of  the  Scattered  Hj>  Number  Density  at.  an  X 
Location  of  0  Meters  For  the  200  Kilometer  Case. 
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Figure  10.  A  Contour  Plot  of  the  Scattered  H-,  Number  Density  at  an  X 
Location  of  500  Meters  For  the  200  Kilometer  Case. 
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Figure  li.  A  Contour  Plot  of  the  Scattered  H2  Number  Density  at  an  X 
Location  of  1500  Meters  For  the  200  Kilometer  Case. 
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Figure  13.  A  Contour  Plot  of  the  Scattered  C02  Number  Density  at  a  2 
Location  of  1000  Meters  For  the  200  Kilometer  Case. 
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NUMBER  DENSITY  (nOLECULES/CMJ) 


’ ,'j.iirft  14.  The  Upstream  Density  Decay  for  the  Three  Species  at  a  Z 
Conation  of  Kilometer  Above  the  Shuttle  for  the  200 
Kilometer  Case. 
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Figure;  15.  A  Contour  Plot  of  the  Z  Velocity  Component  for  the 
Scattered  H2  at  an  X  Location  of  0  Meters  For  the  200 
Kilometer  Case. 
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p’.iyure  16.  A  Contour  F'lot  of  the  Translational  Temperature  for  the 
Scattered  H2  at  an  X  Location  ot  0  Meters  For  the  200 
Kilometer  Case. 
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the  7  km/s  shuttle  velocity),  with  a  minimum  of  approximately  2500  K, 
increasing  to  values  on  the  order  of  6000  K  about  1.5  km  above  the  shuttle. 
The  failure  of  Figures  15  and  16  to  show  precise  riglit-to-lef t  symmetry  is 
aue  to  statistical  scatter,  and  the  degree  of  asymmetry  gives  an  indication 
of  the  error  in  the  calculat ions . 

Figures  17  and  18  show  the  same  two  quantities  for  scattered  C02 .  If 
the  flow  were  in  translational  equilibrium,  of  course,  the  curves  would  be 
identical  for  scattered  Hy  and  scattered  C02 .  It  can  be  seen  that  the 
curves  took  very  different,  indicating  substantial  translational 
nonequilibrium.  The  CO2  reaches  a  maximum  2  veiocity  of  approximately  half 
that  of  the  H2  ,  and  the  maximum  occurs  much  closer  to  the  shuttle.  The 
temperature  for  the  scattered  CU2  is  substantially  larger  than  that  of  H2 , 
however,  reaching  as  high  as  12,000  K.  it  should  be  noted  that  the 
translational  temperature  is  proportional  to  the  molecular  weight  times  the 
square  of  the  velocity  fluctuation  from  the  mean.  Hence,  although  the  C02 
has  a  higher  translational  temperature,  it  corresponds  to  a  lower  velocity 
fluctuation  due  to  the  much  heavier  molecular  weight  for  C02.  The 
translational  temperature,  rather  than  the  velocity  fluctuation,  is 
presented  since  this  is  the  quantity  that  becomes  identical  between  species 
in  equilibrium. 

The  effect  of  molecular  weight  on  the  scattered  species  properties  in 
the  vicinity  of  the  shuttle  (at.  the  origin)  is  illustrated  in  Figures  19  to 
21.  Figure  19  shows  the  normalized  densities  as  a  function  of  molecular 
weight.  The  densities  in  this  figure  have  been  divided  by  the  exhaust  mole 
fraction  for  each  species,  so  the  figure  represents  the  tendency  of  the 
molecules  to  make  it  back  to  the  shuttle  vicinity,  irrespective  of  their 
initial  prevalence.  The  figure  shows  that  the  tendency  to  get  back  to  the 
shuttle  vicinity  is  essentially  independent  of  species  molecular  weight. 
This  somewhat  unexpected  result  Is  apparently  due  to  the  shuttle's 
location.  That  is,  the  outer  reaches  of  the  scattered  cloud  are  relatively 
highly  populated  in  the  more  mobile  H2 .  Since  conservation  applies,  H2  is, 
therefore,  relatively  depleted  in  the  core  region  of  the  plume.  In  between 
the  two  extremes,  there  are  places  where  molecular  weight  does  not  have  a 
great  effect  on  local  density,  arid  the  shuttle  appears  to  be  at  such  a 
place.  The  same  result  was  noted  in  the  calculations  for  the  other 
altitudes  and,  if  this  trend  holds  up,  it  may  facilitate  a  method  of  making 
quick  estimates  of  contar  •  ation  due  to  atmospheric  scattering.  Since  this 
is  probably  the  most  difficult  portion  of  the  contamination  to  calculate, 
it  is  a  potentially  significant  result. 
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Figure  17.  A  Contour  Plot  of  the  L  Velocity  Component  for  the 
Scattered  C02  at.  an  X  Location  of  0  Meters  For  the  200 
Kilometer  Case. 
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Figure  19.  The  Effect  of  Exhaust  Species  Molecular  Weight  on  the 
Normal ized  Scattered  Species  Density  in  the  Vicinity  of 
rhe  Shuttle. 
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Figure  20.  The  Effect  of  Exhaust  Species  Molecular  Weight  on  the 
Scattered  Species  Z  Velocity  Component  in  the  Vicinity  ol 
the  Shuttle. 
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Figure  21.  The  Effect  of  Exhaust  Species  Molecular  Weight  on  the 
Scattered  Species  Translational  Temperature  in  the 
Vicinity  of  the  Shuttle. 
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Figure  20  shows  the  vertical  velocity  component  as  a  function  of 
molecular  weight  at  the  shuttle.  Note  that  the  zero  on  the  scale  is  at  the 
top  of  the  curve,  with  increasingly  negative  values  corresponding  to 
downward  on  the  graph.  Hence,  there  is  a  monotonic  variation,  with  the 
lighter  molecules  showing  a  much  larger  downward  vertical  velocity 
component  than  the  heavier  molecules  when  they  do  make  it  back  to  the 
shuttle.  Finally,  the  translational  temperature  is  shown  in  Figure  21. 
Again,  it  should  be  noted  that  the  temperature  at  the  shuttle  increases  for 
increasing  molecular  weight,  but  that  this  corresponds  to  a  decreasing  mean 
velocity  fluctuation. 

12.3  Surface  Contamination  at  200  Kilometers 

Useful  direct  surface  contamination  statistics  were  not  obtained  for 
this  run.  The  reason  is  that  the  altitude  was  high  enough  that  the 
atmospheric  mean  free  path  was  greater  than  shuttle  dimensions,  ana  the 
solution  region  haa  to  extend  to  many  mean  free  paths  to  describe  the 
contaminant  cloud.  The  resultant  large  cell  size  meant  that  molecular 
interactions  with  the  shuttle  became  sufficiently  improbable  that  none  of 
significance  occurred.  Essentially,  the  problem  results  from  the 
separation  of  the  length  scales  of  the  atmospheric  mean  free  path  ana  the 
shuttle  dimensions.  The  problem  is  even  more  severe,  of  course,  at.  the 
higher  altitudes. 

A  modification  to  the  code  is  planned  whereby  the  present  solution 
wiil  be  calculateo  and  called  the  "outer"  solution.  This  solution  will 
then  define  inward  fluxes  of  all  species  for  an  "inner"  solution  region, 
which  is  just  the  smallest  set  of  the  outer  solution  cells  whii.h  completely 
contain  the  shuttle.  With  the  boundary  condition  given  by  the  outer 
solution,  the  inner  solution  can  then  be  run  with  a  greatly  reduced  grid 
spacing,  and  good  statistics  will  be  obtained  for  direct  surface 
contamination.  The  two  spatial  solutions  can  be  combined  to  give  a 
resultant  cioud  solution  as  well. 

For  the  present  problem,  it  is  possible  to  exploit  the  separation  of 
length  scaies  to  obtain  useful  contaminant  flux  information.  For  a 
Maxwellian  velocity  distribution,  the  one-way  flux,  f + ,  of  molecules  in  any 
direction,  i+,  (that  is,  counting  only  those  molecules  which  have  a 
positive  velocity  component  in  the  specified  direction)  is  given  by: 

f+  =  (n/ (2-ya )  3  { fexp(-W^)  +  W[1  +  erf(W)j)  ,  (136) 
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(137) 


vl  =  -va  U+  (138) 

ana 

L’t  =  u*7  +  (139) 

In  these;  relations,  u  is  the  mean  species  flow  velocity,  m  is  its 
molecular  weight,  n  is  its  number  density,  and  T  is  its  translational 
temperature.  Although  the  actual  species  velocity  distributions  are 
cieariy  non-Kaxwei l ian ,  the  representation  should  be  quite  accurate  for  the 
present,  case  since  the  first  three  moments  of  the  actual  distribution 
function,  with  no  assumption  of  equilibrium,  are  used  to  determine  n,  u  and 
T. 


The  resulting  fluxes  for  the  three  molecules  are  shown  in  Figures  22 
through  24  for  the  Y-Z  plane.  For  scattered  H2 ,  the  flux  reaches  a  maximum 
on  the  top  of  about  1.2xl015  moiecuies/cm2/s ,  which,  if  it  were  to  stick, 
wouia  corresponu  to  approximately  a  monolayer  in  about  half  a  second.  Of 
course,  H2  wouia  not  be  expectea  to  stick,  but  water  would,  ana  the  water 
return  flux  is  70%  as  high  on  the  upper  surface.  The  potential  for 
contamination  via  scacterea  H.^0  would  seem  to  be  substantial  for  this 
example,  ana  CO2  may  also  be  a  problem  although  the  maximum  flux  is 
approximately  a  factor  of  eleven  less  than  for  the  H2O. 

There  is  a  slight  asymmetry  in  the  flux  curves,  which  is  more  apparent 
for  the  scactereo  C02  than  the  others.  This  asymmetry  is  due  to  the 
inclusion  of  the  caicuiateo  Y  velocity  component  in  the  flux  expression. 
From  symmetry,  it  is  known  that  there  would  be  a  zero  Y  velocity  component 
for  ail  species  at  the  origin;  and  any  calculated  value  is,  therefore, 
indicative  of  the  statistical  error.  Although  using  a  zero  Y  velocity 
component  rather  than  the  calculated  one  is  entirely  j  stified  on  physical 
terms,  it  is  useful  to  include  the  calculated  one  precisely  because  it 
gives  a  measure  of  the  error  in  the  calculation. 
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Figure  22. 

The  Return  Flux  of  Scattered  H2  in  the  Y-Z  Plane  as  a 
Function  of  Azimuthal  Angle,  in  the  Vicinity  of  the 

- 

Shuttle,  for  the  200  Km  Case. 
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Figure  23.  The  Return  Flux  or  Scattered  in  the  Y-Z  Plane  as  a 

Function  of  Azimuthal  Angle,  in  the  Vicinity  of  the 
Shuttle,  for  the  200  Km  Case. 


89 
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Figure  24.  The  Return  Flux  of  Scattered  COg  in  the  Y-Z  Plane  as  a 
Function  of  Azimuthal  Angle,  in  the  Vicinity  of  the 
Shuttle,  for  the  200  Km  Case. 
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Another  feature  of  the  flux  curves  which  is  noteworthy  is  that  the 
relative  magnitudes  of  the  maximum  (upper  surface)  and  minimum  (lower 
surface)  fluxes  varies  significantly  with  molecular  weight.  The  heavier 
species  have  a  flux  which  is  more  dominated  by  the  thermal  or  random 
component  of  the  velocity  distribution  function  and  less  so  by  the  mean 
velocity  or  directed  component.  Due  to  this,  the  contaminant  flux  on  the 
lower  surface  for  scattered  COg  is  639s  of  the  upper  surface  flux,  while  for 
scattered  H2  it  is  oniy  27%  of  the  upper  surface  flux. 

12.4  Resuits  at  250  and  300  Kilometers 

The  flux  curves  for  the  250  and  300  km  altitudes  are  given  in  Figures 
25  througn  30.  The  same  general  trends  noted  above  are  evident  at  the 
higher  altitudes  as  well,  except  that  the  level  of  the  fluxes  are  reduced. 
The  most  interesting  resuit  from  the  calculations  at  varying  altitude  is 
the  scaii rtg  of  the  scattered  fluxes.  At  high  enough  altitude,  in  the  free 
molecuiar  flow  regime,  the  return  flux  will  be  due  to  single  collisions 
between  the  exhaust  species  and  an  essentially  undisturbed  atmosphere.  In 
this  flow  regime,  therefore,  the  return  flux  would  be  expected  to  simply  be 
proportional  to  the  ambient  number  density.  As  the  altitude  becomes  lower, 
however,  a  portion  of  the  flux  will  be  due  to  exhaust  species  which  have 
experienced  multiple  collisions  with  atmospheric  species.  The  return  flux 
in  a  multiple  collision  regime  would  be  expected  to  exhibit  a  nonlinear 
behavior  and  show  a  stronger  dependence  on  ambient  number  density. 

Figure  31  shows  the  maximum  returned  fluxes  in  the  Y-Z  plane  as  a 
function  of  ambient  number  density  for  each  of  the  three  species.  Also 
shown  in  the  figure  is  a  dashed  line  showing  the  slope  that  would  result 
from  a  linear  dependence  of  return  flux  on  ambient  number  density.  The 
figure  makes  it  ciear  that  multiple  collisions  are  quite  Important  at  200 
km  (the  highest  ambient,  number  density),  since  the  decline  in  the  maximum 
return  flux  is  much  greater  than  linear  with  ambient  density  when  going  to 
250  km.  The  variation  between  250  and  300  km  is  nearer  the  linear  slope, 
though  still  somewhat  steeper,  indicating  that  multiple  collisions  still 
play  a  role  at  the  250  km  altitude.  One  can  speculate  from  the  figure  that 
if  runs  had  been  made  for  higher  altitudes  (i.e.,  lower  densities?  that  the 
free  molecular  limit  would  be  realized  at  300  km  and  above,  but  this  cannot 
be  proven  in  the  absence  of  such  runs. 
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Figure  25. 

The  Return  Flux  of  Scattered  H2  in  the  Y-Z  Plane  as  a 
Function  of  Azimuthal  Angle,  in  the  Vicinity  of  the 

- 

Shuttle,  for  the  250  Km  Case. 
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MOLECULAR  FLUX  (MOLECULES/CMVS) 
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Figure  26.  The  Return  Flux  of  Scattered  H20  in  the  Y-Z  Plane  as  a 
Function  of  Azimuthal  Angle,  in  the  Vicinity  of  the 
Shuttle,  for  the  250  Km  Case. 
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gure  27.  The  Return  Flux  of  Scattered  COg  in  the  Y-Z  Plane  as  a 
Function  of  Azimuthal  Angle,  in  the  Vicinity  of  the 
Shuttle,  for  the  250  Km  Case. 
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SCATTERED  H2  FLUX  ON  SHUTTLE 
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Figure  28.  The  Return  Flux  of  Scattered  H2  in  the  Y-Z  Plane  as  a 
Function  of  Azimuthal  Angle,  in  the  Vicinity  of  the 
Shuttle,  for  the  300  Km  Case. 
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Figure  29.  The  Return  Flux  of  Scattered  HgO  in  the  Y-Z  Plane  as  a 
Function  of  Azimuthal  Angle,  in  the  Vicinity  of  the 
Shuttle,  for  the  300  Km  Case. 
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Figure  30.  The  Return  Flux  of  Scattered  COg  in  the  Y-Z  Plane  as  a 
Function  of  Azimuthai  Angle,  in  the  Vicinity  of  the 
Shuttle,  for  the  300  Km  Case. 
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AMBIENT  DENSITY  (MOLECULES/CM3) 


Figure  31.  The  Dependence  of  Maximum  Return  Flux  of  the 
Atmospherically  Scattered  Molecules  on  Ambient  Number 
Density  for  the  200,  250  and  300  Km  Runs. 
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13.  CONCLUSIONS 


f.'if  present  version  of  the  SOCRATES  coue  should  be  regarded  as  the 
first.  step  i  n  an  ongoing  development  project.  It  does  not  yet.  have  modules 
tor-  many  of  the  contaminant  sources  discussed  in  the  introduction,  but  it 
was  designed  m  a  modular  fasnion  so  the  inclusion  of  such  additional 
sources  will  be  a  straignt forward  task.  The  mouel  does  describe  what  is 
prooaoiy  cue  most  difficult  portion  of  the  contamination  proulem:  the 
scattering  of  shuttle  piume  molecules  via  the  atmosphere.  it  presents  a 
solid  foundation  for  the  addition  of  new  sources  to  the  model,  which  is 
already  underway. 

The  sample  calculations  presented  here  provide  important  insight  into 
the  rote  of  molecular  weight  on  the  tendency  of  molecules  to  get  scattered 
duck  to  the  shuttle  and  the  level  of  fluxes  to  expect  from  such  scattering, 
it  is  hoped  that  these  calculations,  and  others  to  follow,  will  allow  the 
development  of  simplified  approximate  expressions  which  will  facilitate  the 
estimation  of  contamination  for  general  molecules  and  altitudes.  Such 
estimates  nave  a  substantial  place  in  initial  experiment  design.  More 
complete  calculations,  of  course,  will  always  be  desirable  l'or  experiment 
prediction  ana  data  analysis. 
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APPENDIX  A.  CHI-SQUARE  PROBABILITY  DISTRIBUTION 


A.l  Physical  Basis 

Fundamentally.  the  rhi-square  function  represents  the  distribution  of 
energy  in  an  equilibrium  classical  system  with  v  degrees  of  freedom.  It  is 
a  weii  known  classical  result  that  each  degree  of  freedom  for  a  molecule  in 
an  equilibrium  gas  will  have,  on  the  average,  an  energy  of  kT/2,  where  k  is 
Boltzmann's  constant  and  T  is  temperature.  (For  example,  the  translaiionai 
moue ,  with  three  degrees  of  freedom,  has  an  average  energy  of  3kT/2  per 
molecule.  The  distribution  of  translational  energy  among  the  various 
molecules  follows  a  chi-square  distribution  with  3  degrees  of  freedom.) 
tuner  mooes  of  energy  (molecular  rotation  and  vibration)  have  their  own 
characteristic  number  of  degrees  of  freedom,  which  may  or  may  not  be  fuJ ly 
excneu  in  the  energy  range  of  interest.  If  a  mode  is  not  fully  excited, 
mat  simpiy  means  that  it  is  behaving  as  if  it.  had  a  non-integer  number  ol 
degrees  of  freedom  within  the  classical  approximation.  The  number  of 
internal  degrees  of  freedom  is  directly  related  to  the  heat  capacity  of  the 
gas  and,  essentially,  v  is  selected  to  match  the  known  heat  capacity  of  a 
given  moiecuie  in  a  given  energy  range.  The  assumption  of  a  constant 
numoer  ol  degrees  ot  freedom  is.  therefore,  equivalent  to  the  assumption  of 
a  conscant  heat  capacity.  A  discussion  of  the  implementation  of  such  a 
mooei  allowing  for  a  finite  rate  relaxation  towards  equilibrium  between 
translational  and  internal  mooes  is  given  in  Reference  A-l. 

A. 2  Definition  and  Mathematical  PropertiesA~2 

The  chi-square  probability  density  function,  f(X;v),  defines  a 
distribution  ol  X  in  a  domain  of  zero  to  infinity  via 


f (X;v) 


x(v/2  -  1 >cxp(-X/2 ) 
2  (  V/2  )  p(  -y/2  ) 


(A-l) 


A-:.  Elgin,  J.  B .  ,  "Getting  the  Good  Bounce:  Techniques  for  Efficient 
Monte  Cano  Analysis  ot  Complex  Reacting  Flows,"  Report  SS1-TR-28, 
Spectral  Sciences,  inc.,  Burlington,  MA  (1983). 

A-2.  Aoramowitz.  M.  ano  Stegun,  i.  A.,  Handbook  of  Mathematical 
Functions,  National  Bureau  ot  Standards,  940  (1968). 


where  v  is  a  positive  parameter  of  the  distribution  referred  to  as  the 
number  of  decrees  of  freedom.  The  chi-square  distribution  results  in  a 
mean  value  of  X  equal  to  y.  Figure  A-l  is  a  plot  of  the  chi-square 
probability  density  function  for  v  equal  to  1,  2,  and  3. 

The  chi-square  distribution  has  a  fundamental  addition  property  such 
chat  if  Xj  is  selected  from  a  chi-square  distribution  with  v,  degrees  of 
freedom,  ana  X2  is  selected  from  a  chi-square  distribution  with  v 2  degrees 
of  freedom,  then  their  sum  will  be  distributed  according  to  a  chi-square 
distribution  with  v.  +  v.,  degrees  of  freedom.  This  property  is  of 
substantial  theoretical  arid  practical  importance. 

Ff  the  variable  Z  is  distributed  according  to  a  normal  distribution 
with  zero  mean  and  unit  variance,  then  will  follow  a  chi-square 
distribution  with  one  uegree  of  freedom.  Tt  follows  from  the  above 
addition  property  that,  in  general  ,  if  Zj,  Z2,  •••  ,  Zn  are  n  variables 
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Figure  A-l.  A  Plot  of  the  Chi-Square  Probability  Density  Function 
for  v  Equal  to  1,  2,  arid  3. 
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seiecten  from  a  such  a  normal  distribution,  and  X  is  defined  as  the  sum  of 
the  squares  of  the  Z j,  then  the  X's  that  result  will  be  distributed 
according  to  a  chi-square  distribution  with  n  degrees  of  freedom. 

Finally,  if  t  is  distributed  according  to  a  probability  density 
func  i.ion  g(t;p,q),  where 


« ( t ;  p ,  q ) 


tP-l(1  -  tj^-1 
B(P,q) 


( A-2 ) 


and 


H  ( p .  q ) 


1 

=  |  tP-1(l  -  t)<l_1dt. 

J 

0 


npjnq] 

r(p  +  q) 


(B  is  the  Beta  function)  then  t  can  be  sampled  via 


( A-3 ) 


t 


Xi  ♦ 


(A-4) 


where  X1  is  selected  from  a  chi-square  distribution  with  Vj  degrees  of 
freedom,  and  X2  is  selected  from  a  chi-square  distribution  with  v2  degrees 
of  freedom,  with 

vl  =  2p  .  (A-5) 


and 


v2  =  2q  ( A-6 ) 

The  significance  of  Eq.  (A-4)  is  that  it  reduces  the  sampling  from  a  two- 
parameter  distribution  (F.q.  (A-2))  to  two  samplings  from  a  one-parameter 
distribution.  The  distribution  represented  by  Eq.  (A-2)  arises  in  cases 
wriere  a  constrained  amount  of  total  energy  is  distributed  among  various 
mooes,  and  its  relation  to  the  chi-square  distribution  apparently  has  not 
been  appreciated  by  developers  of  techniques  for  Monte  Carlo  fluid 
mccnanics . 
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A. 3  Sampling  From  a  Chi-Square  Distribution 

The  neea  for  sampling  from  a  chi-square  distribution  comes  up  when 
sampling  initial  values  of  internal  energies,  when  calculating  inelastic 
collisions  via  the  statistical  collision  modelA"^  or  when  calculating  the 
equilibrium  aftermath  of  many  collisions  in  a  cell.  Since  these  operations 
must  be  performed  repeatedly  in  the  heart  of  a  Monte  Carlo  simulation,  it 
is  important  that  the  sampling  be  done  efficiently  and  accurately. 

For  clarity,  the  result  of  each  sampling  method  discussed  beiow  will 
oe  denoted  by  a  different  letter  subscript  to  X.  All  sampling  procedures 
make  use  of  a  random  number  generator  which  returns  a  number,  R,  selected 
from  a  probability  density  which  is  uniform  on  the  interval  between  zero 
ano  one.  Each  occurrence  of  R  indicates  a  distinct  sampling  from  the 
random  number  generator. 

A. 3.1  Analytic  Sampling  for  Integer  v 

Direct  sampling  of  Eq.  (A-l )  can  be  performed  for  integer  v,  as  shown 
below. 


A. 3. 1.  1  V  —  0 

As  u  (an  intrinsically  non-negative  quantity)  approaches  zero,  the 
distribution  function  approaches  a  delta  function,  and  a  proper  sampling  is 
achieved  by  simply  selecting 

=  0  (A-7) 

A. 3.  i  ,‘i  v  =  1 

For  sampling  with  v  =  1  (as  well  as  for  several  other  cases),  it  is 
convenient  to  introduce  the  transformation  Z^  =  X.  Z  is  then  distributed 
according  to  the  probability  density  function  p ( Z )  given  by 


A-3.  Borgnakke,  C.  and  Larsen,  P.  S.  ,  "Statistical  Collision  Model  for 
Monte  Carlo  Simulation  of  Polyatomic  Gas  Mixture,"  Journal  of 
Computational  Physics,  18 ,  405  (1975). 
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Z( V-1 jexp(-Z2/2) 
(V/2  -  l)p(v/2) 


(A-8) 


For  v  =  1,  this  distribution  is  simply  a  normal  distribution  adjusted  to 
allow  for  positive  only  argument.  Sampling  from  this  distribution  is 
described  in  Reference  A-4 .  When  the  result  is  cast  back  in  terms  of  X, 
the  result  is 


A  =  2trR 


(A-9 ) 


ana 


X„  =  -2log(R)sin2(A)  (A-10) 

A. 3. 1 .3  v  =  2 

Wnen  v  =  2,  the  integral  of  Eq.  (A— 1 )  can  be  analytically  inverted, 
leading  to  the  direct  sampling 

Xc  =  -2iog(R)  (A- 11 ) 

A. 3.1. 4  v  Equal  to  an  Even  Integer 

The  extreme  simplicity  of  the  above  sampling  for  v  =  2,  together  with 
the  addition  property  of  the  chi-square  distribution,  means  that  sampling 
tor  v  equal  to  an  even  integer  is  quite  direct.  Let  J  *=  v/2,  then  a  proper 
chi-square  sampling  is  given  by 

Xa  =  -210g(R1R2. . .Rj)  ,  (A-12 ) 

where  Rj  througn  Rj  denote  J  samplings  from  the  random  number  generator. 
The  fact  that  the  log  need  only  be  taken  once  in  Eq.  (A-12)  means  that  the 
evaluation  of  XQ  is  quite  efficient,  even  for  moderately  large  v. 


A-l  Bird,  G.  A.,  Molecular  Gas  Dynamics.  Clarendon  Press,  Oxford 
(1976). 
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A. 3. 3.5  v  Equal  To  an  ood  integer 


For  v  equal  to  an  odd  integer,  the  addition  property  of  the  chi-square 
distribution  allows  the  simple  combination  of  the  results  for  v  equal  to 
one  anu  v  equal  to  an  even  integer,  i.e., 

xe  =  xb  T  xd  ■  ( A— 13) 

where  Xh  is  given  in  Eq .  (A-10)  ana  Xd  is  given  in  Eq.  (A-12)  with  J  = 
(V-1J/2. 

A.  3. 2  Generalized  Acceptance-Rejection  Sampling 

For  non-integer  v,  it  is  necessary  to  use  a  general  izea  form  of 
acceptance-rejection  sampling.  Before  the  application  to  chi-square 
sampling  is  presented,  the  acceptance-rejection  technique  arm  its 
generalization  will  be  briefly  discussed. 

A. 3. 2.1  Starmara  Acceptance-Rejection  Sampling 

The  usual  acceptance-rejection  technique  for  sampling  from  a  general 
distribution  function,  p(x),  proceeds  as  follows: 

1)  The  domain  of  x  is  approximated,  if  necessary,  by  a  finite 
sub-domain. 

£ 

2)  The  maximum  value  of  p(x) ,  p  ,  is  calculated. 

3)  A  variable  %  is  selected  from  the  domain  of  x  via 

^  =  xmin  +  R(xmax  xmin^ 

4)  p(£)/p*  is  calculated,  and  another  random  variable,  R,  is 
generated  x  is  set  equal  to  f;  if  R  is  less  than  p(^)/p  . 

5)  Steps  3  and  4  are  repeated  until  a  value  of  x  is  determined. 

Note  that  the  probability  of  acceptance  of  the  random  variable  in  step  4  is 
proportional  to  the  distribution  function  being  sampled,  so  the  resulting  x 
values  will  follow  the  desired  distribution  function. 
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Although  the  generality  of  this  approach  makes  it  very  powerful,  it. 
does  suffer  from  the  following  drawbacks: 

•  If  the  distribution  function  differs  significantly  from  its 

maximum  value  within  a  substantial  portion  of  the  sampled  domain, 
then  the  rejection  rate  may  be  high-  This  obviously  leads  to  a 
siow  sampling  procedure. 

•  If  the  finite  sub-domain  is  reduced  to  increase  the  acceptance 

rate,  then  the  sampling  deviates  from  the  true  distribution 
function. 

•  The  procedure  is  incapable  of  sampling  from  an  unbounded 

distribution  function. 

A. 3. 2. 2  Generalization  of  the  Acceptance-Rejection  Technique 

The  following  procedure  comprises  a  generalization  of  the  acceptance- 
rejection  technique: 

1)  A  second  distribution  function,  q(x),  which  can  be  sampled 

analytically  is  chosen.  Conuitions  on  q(x)  will  be  discussed 
be  low. 

2)  The  maximum  value  of  p(x)/q(x),  (p/q)*,  is  calculated. 

3)  A  variable,  £,  is  sampled  from  q. 

£ 

4 )  lp(S)/q(S) j/(p/q)  is  calculated,  and  another  random 
variable,  R,  is  generated,  x  is  set  equai  to  5  if  R  is  less  than 
Q 

5)  Steps  3  ana  4  are  repeated  until  a  value  of  x  is  determined. 

It.  snouid  be  noted  that  the  probability  density  for  a  given  value  of  x 
is  proportional  to  the  product  of  the  initial  selection  probability  times 
tne  acceptance  probability.  Since  the  former  probability  is  proportional 
to  q ( x ) ,  ana  tne  latter  is  proportional  to  p(x)/q(x),  the  distribution  of 
accepted  values  does  indeed  fol low  the  distribution  function  p(x). 

Fhe  usual  acceptance-rejection  technique  is  simply  the  case  where  q(x) 
is  constant,  but  it  is  evident  that  this  is  not  always  the  best  (or  even  a 
possioie)  choice.  All  of  the  objections  to  the  standard  acceptance- 
rejection  tecnnique  can  be  removed  or  ameliorated  by  a  suitable  choice  for 
q  1  x  )  Iri  particular: 


•  There  is  no  neeu  to  approximate  the  domain  of  x  with  a  finite 
sub-domain.  ft  is  merely  necessary  that  the  domain  for  q  include 
the  domain  for  p.  The  domain  for  q  can  be  larger  than  that  for  p, 
since  whenever  a  value  is  selected  from  outside  the  domain  for  p 
it  will  always  be  rejected  in  step  4  above. 

•  If  q  is  selected  to  bo  close  to  p,  at  least  in  the  region  of 
highest  probability,  then  the  acceptance  rate  of  trial  values  will 
be  large. 

•  Unbounded  distribution  functions  can  be  sampled  if  q  is  chosen  to 
have  the  same  type  of  singularity  as  p,  since  the  only  requirement 
is  that  the  ratio  (p/q)  remain  bounded. 

For  any  given  situation,  the  choice  of  the  function  q  is  a  bit  of  an 
art,  guided  by  the  concerns  highlighted  above:  q  must  have  a  domain  which 
includes  the  domain  of  p;  p/q  must  remain  bounded;  and  (p/q)  should  achieve 
its  maximum  in  the  vicinity  of  the  maximum  of  p. 

A. 3. 3  Exact  Acceptance-Rejection  Sampling  for  a 
Chi-Square  Distribution  with  Large  v 

The  acceptance-rejection  technique  described  above  can  be  used  to 

achieve  an  exact  sampling  from  a  chi-square  distribution  for  large  v. 

(Actual ly,  the  approach  is  perfectly  valid  for  all  v  >  1,  but  the  method  to 

be  described  in  Subsection  A. 3. 5  is  to  be  preferred  for  v  <  45,  or  so.) 

The  procedure  utilizes  the  transformed  chi-square  distribution,  p ( 2 ) ,  given 

by  Eq.  (A-8)  as  the  distribution  to  be  sampled.  A  normal  distribution  is 

used  as  the  initial  distribution  which  can  be  sampled  analytically.  The 

normal  distribution  is  chosen  to  have  a  unit  variance  and  a  mean  which 

corresponds  to  the  location  of  the  maximum  of  p(Z).  This  maximum  occurs  at 
* 

Z  given  by 

Z*  =  -yv  -  1  .  (A-14 ) 

The  functional  form  of  the  normal  distribution,  q(Z),  is 


q(Z)  =  exp[-|(Z  -  Z*)2]/t/2it 


( A-15 ) 
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union  not  only  Has  a  maximum  at  the  same  .location  as  Eq.  (A-8),  but  has  the 
same  exponential  factor  as  2  approaches  infinity  arul  a  domain  which 
mu maos  tn at  of  Eq.  (A-8).  The  sampling  of  a  chi-square  value  proceeds  as 
fo i tows : 

'  **■  * —  — — 

1)  2  =  tv  -  l  is  calculated. 

2)  A  sample  from  the  distribution  given  by  Eq.  (A-15)  is  taken  via: 


a)  A  =  2trR  (A-16) 

b)  8  =  -Jog(R)  (A-17) 

c)  Z  =  l*  +  -*~2B  sin  (A)  (A-18) 

3)  The  acceptance  probability,  Q,  is  computed  as 

Q  =  (Z/Z*)(v  “  1,exp[-Z*(Z  -  Z*)]  .  (A-19) 

(0  is  taken  to  be  zero  for  negative  Z.) 


4)  Another  random  variable  is  generated,  and  Z  is  kept  if  Q  >  R.  If  Z 
is  rejected,  then  steps  2-4  are  repeated  until  a  Z  value  is 
accepted . 

5)  When  a  Z  value  is  accepted,  then  the  corresponding  chi-square 
value  is  given  by 

Xe  =  Z2  (A-20) 

This  procedure  is  illustrated  in  Figure  A-2  which  shows  p(Z),  q(Z)  and  Q ( Z ) 
for  v  =  50 .  Note  that  the  acceptance  probability  is  near  unity  in  the 
vicinity  or  the  maxima  of  the  two  distribution  functions,  so  a  large 
traction  of  the  selected  samples  of  q  ( Z )  will  be  accepted  as  samples  of 
PlZ). 


A. 3.4  Exact  Acceptance-Rejection  Sampling  for  a  Chi-Square 
Distribution  with  (0  <  v  <  2). 

For  this  domain  of  v  it  is  convenient  to  introduce  another 
transformation  to  Eq.  ( A- 1 ) .  If  W  is  defined  by 
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FUNCTIONAL  VALU 


C\J 


3.0  5.0  7.0  9.0  11.0 


SAMPLED  VARIABLE.  Z  =  X'5 


Figure  A-2 .  A  Representation  of  the  Transformed  Chi-Square 
Distribution,  p(Z),  for  v  =  50  (solid  line).  p ( Z )  is 
sampled  by  first  selecting  a  variable  from  the  shifted 
normal  distribution  (dashed  line)  and  keeping  it  with  a 
probability  given  by  Q(Z)  (dotted  line). 
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exp( -X/2 ) 


(A-21  ) 


then  the  probability  density  function  for  W  is  given  by  h(W),  where 


n  ( W ) 


f og(W) 


( v/2  -  1  ) 


r(v/2! 


( A-22 ) 


The  domain  for  W  is  finite  (between  0  and  1),  but  h(W)  becomes  infinite  us 
W  approaches  unity.  Tito  generalized  acceptance-rejection  technique  can 
stUi  be  used,  however,  since  the  function  q(W)  given  by 

q(W)  =  ( v / 2 ) ( 1  -  W)(v/2  "  1}  (A-23) 

nas  the  same  type  or  singularity  and  can  be  analytically  sampleo.  The 
chi-square  sampling  for  (0  <  v  <  2)  proceeds  as  follows: 

1}  A  sample  from  q(W)  is  generated  via 

W  =  1  -  (A-24) 

2)  The  acceptance  oronability,  Q,  is  computed  from 

Q  =•  l (W  -  1)/ log(W) ] (1  *  v/2)  (A-25) 

3)  Another  random  variable,  R,  is  generated,  and  W  is  kept  if  Q  >  R: 
otherwise  steps  2  and  3  are  repeated  until  a  value  for  W  is 
accepted . 

4)  When  a  vaiue  for  W  is  accepted,  the  corresponding  chi-square  value 
is  given  by 

Xf  =  -2!og(W)  (A-26) 

This  procedure  is  illustrated  in  Figure  A-3,  which  shows  the  two 
distribution  functions,  h(W)  and  q(W)  ,  and  the  acceptance  probability. 
U(^).  tor  v  =  1.  [t  can  be  seen  that  q(W)  provides  an  excellent  choice  for 
the  initial  selection  of  W,  since  the  acceptance  probability  remains  high 
throughout  tne  important  domain  of  W.  This  point  will  be  discussed  in  more 
u« ta  l  I  in  See .  A .  3  .  b  . 
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FUNCTIONAL  VALUE 


SAMPLED  VARIABLE.  W=EXP(-X/2) 


Sigure  A-3 .  A  Representation  of  the  Transformed  Chi-Square 
Distribution,  h(W),  for  v  =  1  (solid  line).  h(W)  is 
sampled  by  analytically  selecting  a  variable  from  q(W) 
(dashed  line)  and  keeping  it  with  a  probability  given 
by  Q(W)  (dotted  line). 
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A. 3. 5  Exact  Chi-Square  Sampling  for  General  v 


Using  the  fundamental  audition  property  of  chi-square  distributions, 
it  is  possible  to  combine  the  proceuure  described  in  Subsection  A. 3. 1.4  for 
v  equal  to  an  even  integer  with  the  procedure  described  in  Subsection  A. 3.4 
tor  ( 0  <  v  <  2 )  to  achieve  a  exact  general  sampling  technique  for  arbitrary 
v.  inis  is  given  simply  by 

X,,  -  X„  *  Xf  ,  (A-27) 

where  Xu  is  calculated  from  Eq.  (A--12)  with  J  equal  to  the  integer  portion 
oi  y/2,  ano  Xf  is  calculated  as  in  the  preceding  section  with  v  being 
replaced  oy  v-2 J . 

It  is  to  be  noted  that  both  the  approach  given  in  this  subsection 
(Eq.  (27))  and  that  given  in  Subsection  A. 3. 3  (Eq.  (A-20))  are  exact  and 
applicable  for  v  >  i.  In  general,  the  approach  of  this  subsection  is 
considerably  faster,  although  as  v  gets  large  the  approach  of  Subsection 
A. 3. 3  becomes  more  attractive.  There  are  two  potential  difficulties  with 
Eq.  ( A - 2  7 )  as  v  becomes  very  large.  Firstly,  the  product  required  in 
Eq.  (A-12)  gets  more  ano  more  cumbersome  to  compute  as  v  increases  and, 
secondly,  the  larger  the  number  of  factors  in  this  product  the  greater  is 
me  chance  that  it  will  yield  a  number  so  small  as  to  produce  a  floating 
point  underflow  on  a  computer.  (Since  Monte  Carlo  codes  must  be  highly 
reliable,  any  such  problem  should  be  made  essentially  impossible.)  It 
turns  out  that  the  second  problem  is  more  restrictive  (at  least  for  32-bit. 
computers),  dictating  that,  the  Eq.  (20)  should  be  used  for  v  greater  than 
45  or  so.  This  Keeps  the  probability  of  an  underflow  below  10~10  on  any 
given  sampling. 

A. 3.6  Approximate  Chi-Square  Sampling  for  (0  <  v  <  2) 

Trie  procedure  described  in  the  preceding  section  is  quite  efficient, 
mii.  it  is  nonetheless  useful  to  consider  approximate  methods  for  sampling 
from  crti-square  distributions.  While  it  would  be  scarcely  possible  to 
improve  on  sampi ing  for  even  integer  v  discussed  m  Subsection  A. 3. 1.4,  it 
is  ’-easonable  to  investigate  approximations  for  the  (0  <  v  <  2)  portion 
discussed  in  Subsection  A. 3. 4.  A  likely  piace  to  look  for  useful 
approximations  in  this  proceuure  is  m  the  calculation  of  Q  (Eq.  (A-25)), 
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wmen  must  be  performed  for  every  W  selected  in  Eq.  (A-24).  (Noie  that  the 
cai  cm  aliort  or  (J  involves  more  computational  effort  than  the  calculation  of 
W.  ) 

The  overall  probability  that  the  value  chosen  in  Eq.  (A-24)  wi  i  i  be 
xept  as  a  sample  of  Eq.  (A-22)  Is  given  by  P,  where 


r 

p  =  |  q ( w ) g ( w ) aw  =  r( 1  +  v/2) 

J 

o 


(A -20) 


Hence,  as  v  approaches  0  or  2,  all  initially  selected  values  of  W  are  kept 
as  valid  samples  of  Eq.  { A— 22 } .  and  the  computation  of  Q  serves  no  useful 
purpose.  In  the  worst  case  (v  =  0.92)  the  overall  acceptance  probability 
is  89%,  ana  only  11%  of  the  initially  selected  variables  are  rejected.  The 
approximate  chi-square  sampling  involves  approximating  Q(W)  by  an  easily 
calculable  function  which  differs  little  from  Eq.  (A-25).  The  current 
approximation  is  Qa(W),  given  by 


Qa(W)  =  1  -  (1  -  V/2 ) ( 1  -  W) ( . 5  +  a( v) ( 1  -  W)2j  ,  (A-29) 


wnere 


a( V )  =  .2511V  +  .2073  (A-30) 

Qa  was  selected  to  match  the  value  ana  slope  of  Q  at  W  =  1 ,  which  is  the 
region  of  nighest  probability  uensity.  The  coefficient  <x(v)  is  a  linear 
tit  to  values  chosen  to  be  optimal  in  the  least  squares  sense.  A 
comparison  of  Q  arid  Qa  for  v  =  1.0  is  shown  in  Figure  A-4 ,  which 
demonstrates  the  suostantiai  accuracy  of  the  approximation. 

It  is  f uiiaamentaily  more  important,  ot  course,  to  compare  the  correct 
cm -square  attribution  with  the  distribution  which  is  effectively  being 
sampled  m  the  approximate  technique.  If  h3(W)  is  the  approximation  analog 
ot  h ( W ) ,  then  ha(W)  is  proportional  to  the  product  q(W)Q(W),  i.e., 

ha(W)  =  A(l-W)(v/2  "  ^Q^W)  .  (A-31) 

where  the  normalization  factor,  A,  is  determined  by  requiring  that  ha(W) 
give  unity  when  integrated  between  zero  and  one.  This  results  in 
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v2  r  2v  t  8  2a( V ) 

2v^  +  4V  v  -»  6 


(A-32) 


Once  ha(W)  is  defined,  the  corresponding  distribution  function  for  X, 
fa(X;v)  is  obtained  by  multiplying  ha(W)  by  the  magnituue  of  dW/dX  (=  W/2), 
anu  substituting  W  =  exp(-X/2).  The  comparison  between  f(X;1.0)  and 
f  a  ( X ;  1 . 0 )  is  given  in  Figure  A-5,  ano  the  agreement,  is  excellent.  The  use 
of  the  approximate  technique  is  approximately  40%  faster  than  the  exact 
acceptance-rejection  technique,  and  the  difference  in  the  distributions 
oemg  sampled  will  probably  always  be  negligible.  Although  the  ability  to 
sample  from  an  exact  chi-square  distribution  will  be  kept  as  an  option,  it 
is  felt  that  the  approximate  technique  offers  a  substantial  time  savings 
for  an  inconsequent lal  loss  of  accuracy. 
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ACCEPTANCE  PROBABILITY 


Figure 
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A-4 .  A  Comparison  or'  the  Exam,  and  Approximate  Acceptance 
Probabilities  for  v*l . 
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CHI  SQUARE  PROBABILITY  DENSITY 


Figure 
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A -5.  A  Comparison  ol'  the  Exact  Chi-Square  Distribution, 
t(X;v)  with  the  Approximate  Distribution,  i’a(X;v)  which 
is  Effectively  Being  Sampled  by  the  Approximate 
Technique  Presented  in  Subsection  A. 3. 6.  v  was  taken 
to  be  unity  for  the  comparison. 


A-17 


A. 4  APPENDIX  A  REFERENCES 

A-l.  Elgin,  J.  B.  ,  “Getting  the  Good  Bounce:  Techniques  for  Efficient 
Monte  Carlo  Analysis  of  Complex  Reacting  Flows,"  Report  SSi-TR-28,  Spectral 
Sciences,  Inc.,  Burlington,  MA  (1983). 

A-2.  Abramowitz,  M.  and  Stegun,  I.  A.,  Handbook  of  Mathematical  Functions. 
National  Bureau  of  Standards,  940  (1968). 

A-3.  Borgnakke,  C.  and  Larsen,  P.  S.,  "Statistical  Collision  Model  for 
Monte  Carlo  Simulation  of  Polyatomic  Gas  Mixture,*'  Journal  of  Computational 
Physics,  18,  405  (1975). 

A-4 .  Bird,  G.  A.,  Molecular  Gas  Dynamics,  Clarendon  Press,  Oxford  (1976). 


A-18 


